function load_GHCN_station_list, restore=restore

if keyword_set(restore) then begin
    restore,file='v2.temperature.inv.txt.sav'
    return, s
endif

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
name = STRCOMPRESS(name)
s = { GHCN_station_list, $
      icc: '', wmo: '', wmo_mod: '', name: '', $
      lat: 0.0, lon: 0.0, $
      ele: 0, tele:0, p:'', pop:0, tp:'', v:'', lo:'', co:0, $
      a:'', ds:0, veg:'', bi:'' }
s = REPLICATE( s, n_elements(icc) )
s.icc = icc & s.wmo = wmo & s.wmo_mod = wmo_mod & s.name=name
s.lat = lat & s.lon = lon
s.ele = ele & s.tele = tele & s.p = p & s.pop=pop & s.tp=tp & s.v=v
s.lo = lo & s.co=co & s.a=a & s.ds=ds & s.veg=veg & s.bi=bi
save, s,file='v2.temperature.inv.txt.sav'
return, s
end


function load_GISS_station_list, restore=restore

if keyword_set(restore) then begin
    restore, file='station_list.txt.sav'
    return, s
endif

file='station_list.txt'
fmt = '          A5, A3,     A1, A1,     A32, A35,  A3'
readfmt,file,fmt,wmo,wmo_mod,n,  dummy,  name,dummy,cc, skipline=1
name = STRCOMPRESS(name)
wmo = STRING(wmo,FORMAT='(I5.5)')
wmo_mod = STRING(wmo_mod,FORMAT='(I3.3)')
cc = STRING(cc,FORMAT='(I3.3)')
s = { GISS_station_list, wmo:'', wmo_mod: '', n:'', name:'', icc:'' }
s = REPLICATE( s, n_elements(n) )
s.wmo = wmo & s.wmo_mod=wmo_mod & s.n=n & s.name=name & s.icc=cc
save, s, file='station_list.txt.sav'
return, s
end


function load_station_data, giss, ghcn, restore=restore

if keyword_set(restore) then begin
    restore, 'station_data.sav'
    return, s
endif

ns = n_elements(GISS)
nyr = 2010-1880
arr = fltarr( nyr )
s0 = { months, data:fltarr( nyr, 12 ), $
       names:['January','February','March','April','May','Jun',$
              'July','August','September','October','November','December'] }
s1 = { seasons, data:fltarr( nyr, 4 ), $
       names:['DJF','MAM','JJA','SON'] }
s = { station_data, $
      icc: '', wmo: '', wmo_mod: '', name: '', $
      years: arr, $
      annual: arr, $
      months: s0, $
      seasons: s1 }
s = replicate( s, ns )

for i = 0, ns-1 do begin
;for i = 5308, 5319 do begin ; testing
;for i = 5308, ns-1 do begin ; testing
    station = giss[i].icc + giss[i].wmo + giss[i].wmo_mod + giss[i].n
    file = FINDFILE('./data/'+station+'*.txt')
    if n_elements(file) ne 1 then stop else file = file[0]
    if file eq '' then continue
    readcol, file, year,jan,feb,mar,apr,may,jun,jul,aug,sep,oct,nov,dec,djf,mam,jja,son,ann, skipline=1, /silent
    s[i].icc = giss[i].icc
    s[i].wmo = giss[i].wmo
    s[i].wmo_mod = giss[i].wmo_mod
    s[i].name = giss[i].name

    s[i].years = year
    s[i].months.data[0:n_elements(year)-1,0:11] = $
      [[jan],[feb],[mar],[apr],[may],[jun], $
       [jul],[aug],[sep],[oct],[nov],[dec]]
    s[i].seasons.data[0:n_elements(year)-1,0:3] = [[djf],[mam],[jja],[son]]
    s[i].annual = ann
    ygd = where( s[i].years ne 0, count )
    if count ne 0 then begin
       y0 = min( s[i].years[ygd] )
       s[i].years = SHIFT( s[i].years, y0-1880 )
       s[i].months.data = SHIFT( s[i].months.data, y0-1880 )
       s[i].seasons.data = SHIFT( s[i].seasons.data, y0-1880 )
       s[i].annual = SHIFT( s[i].annual, y0-1880 )
    endif

    print, i
endfor
save, s, filename='station_data.sav'
return, s
end







pro load_panoply,r,g,b
openr, lun, 'panoply_diff.pa1', /get
r = ( g = ( b = bytarr(256) ) )
a = assoc(lun, r)

; load to IDL
r = a[0] & g = a[1] & b = a[2]
r[0] = ( g[0] = ( b[0] = 0 ) )
r[255] = ( g[255] = ( b[255] = 255 ) )
r[254] = ( g[254] = ( b[254] = 192 ) )
tvlct,r,g,b

; return basic color vectors, 8 blue 8 red
r = a[0]
g = a[1]
b = a[2]
subset = makex(8,256,16)
r = r[subset]
g = g[subset]
b = b[subset]
free_lun, lun
end





function make_meta_data, giss, ghcn, data, restore=restore
if keyword_set(restore) then begin
    restore, 'meta.sav'
    return, s
end
s = { wmo: '', wmo_mod: '', icc: '', $
      baseline: 0.0, nbaselineyr: 0, $
      mall: 0.0, ball: 0.0, nall: 0, $
      m100: 0.0, b100: 0.0, n100: 0, $
      m50: 0.0, b50: 0.0, n50: 0, $
      m25: 0.0, b25: 0.0, n25: 0, $
      m2000: 0.0, b2000: 0.0, n2000: 0, $
      m1990: 0.0, b1990: 0.0, n1990: 0, $
      m1980: 0.0, b1980: 0.0, n1980: 0, $
      m1970: 0.0, b1970: 0.0, n1970: 0, $
      m1960: 0.0, b1960: 0.0, n1960: 0, $
      m1950: 0.0, b1950: 0.0, n1950: 0 }
s = replicate(s,n_elements(giss))
for i=0, n_elements( giss )-1 do begin

    giss_idx = i & print, i
    s[i].wmo = giss[giss_idx].wmo
    s[i].wmo_mod = giss[giss_idx].wmo_mod
    s[i].icc = giss[giss_idx].icc

    ghcn_idx = where( giss[giss_idx].wmo eq ghcn.wmo AND $
                      giss[giss_idx].wmo_mod eq ghcn.wmo_mod AND $
                      giss[giss_idx].icc eq ghcn.icc, n )
    if n ne 1 then stop else ghcn_idx = ghcn_idx[0]
    data_idx = where( giss[giss_idx].wmo eq data.wmo AND $
                      giss[giss_idx].wmo_mod eq data.wmo_mod AND $
                      giss[giss_idx].icc eq data.icc, n )
    ;if n ne 1 then stop else data_idx = data_idx[0]
    if n eq 0 then continue
    if n ne 1 then data_idx = data_idx[0]

    ; 1951 to 1980 baseline
    baseyears = where( data[data_idx].years ge 1951 and $
                       data[data_idx].years le 1980 AND $
                       data[data_idx].years ne 0 AND $
                       data[data_idx].annual ne 999.9, n )
    if n ge 2 then begin
        s[i].baseline = total(data[data_idx].annual[baseyears]) / float(n)
        s[i].nbaselineyr = n
    endif

    ; 1st order line fit for all available years
    gd = where( data[data_idx].years ne 0 AND $
                data[data_idx].annual ne 999.9, n )
    if n ge 2 then begin
        s[i].nall = n
        line = indgen(n_elements(data[data_idx].years))
        r1 = poly_fit( line[gd], data[data_idx].annual[gd], 1 )
        s[i].ball = r1[0]
        s[i].mall = r1[1]
    endif

    ; 1st order line fit for last 100 years
    gd = where( data[data_idx].years ne 0 AND $
                data[data_idx].annual ne 999.9 AND $
                data[data_idx].years GE 1908, n )
    if n ge 2 then begin
        s[i].n100 = n
        line = indgen(n_elements(data[data_idx].years))
        r1 = poly_fit( line[gd], data[data_idx].annual[gd], 1 )
        s[i].b100 = r1[0]
        s[i].m100 = r1[1]
    endif

    ; 1st order line fit for last 50 years
    gd = where( data[data_idx].years ne 0 AND $
                data[data_idx].annual ne 999.9 AND $
                data[data_idx].years GE 1958, n )
    if n ge 2 then begin
        s[i].n50 = n
        line = indgen(n_elements(data[data_idx].years))
        r1 = poly_fit( line[gd], data[data_idx].annual[gd], 1 )
        s[i].b50 = r1[0]
        s[i].m50 = r1[1]
    endif

    ; 1st order line fit for last 25 years
    gd = where( data[data_idx].years ne 0 AND $
                data[data_idx].annual ne 999.9 AND $
                data[data_idx].years GE 1983, n )
    if n ge 2 then begin
        s[i].n25 = n 
        line = indgen(n_elements(data[data_idx].years))
        r1 = poly_fit( line[gd], data[data_idx].annual[gd], 1 )
        s[i].b25 = r1[0]
        s[i].m25 = r1[1]
    endif


    ; every decade!
    for yr=1950,2000,10 do begin
       gd = where( data[data_idx].years ne 0 AND $
                   data[data_idx].annual ne 999.9 AND $
                   data[data_idx].years GE yr AND $
                   data[data_idx].years LE yr+9, n )
       if n ge 2 then begin
          str = 's[i].n'+STRTRIM(yr,2)+'=n' & dummy=execute(str)
          line = indgen(n_elements(data[data_idx].years))
          r1 = poly_fit( line[gd], data[data_idx].annual[gd], 1 )
          str = 's[i].b'+STRTRIM(yr,2)+'= r1[0]' & dummy=execute(str)
          str = 's[i].m'+STRTRIM(yr,2)+' = r1[1]' & dummy=execute(str)
       endif
    endfor
 endfor

save, s, filename='meta.sav'
end






pro make_img, s
; s stands for 'station'

sname = s.icc + s.wmo + s.wmo_mod
if sname eq '' then return

set_plot,'ps'
device, /color, bits=8, filename='images/'+sname+'.ps'
load_panoply

nyrs = n_elements(where( s.years ne 0 ) )
if nyrs le 20 then MESSAGE, "Less than 20 data points at this station", /CONTINUE

dd = s.months.data
bad = where( dd eq 0 or dd eq 999.90, nbad, complement=gd, ncomplement=ngd )
if n_elements(gd) le 2 then begin
   MESSAGE, "Error... Not enough data to produce image", /CONTINUE
   return
end
img = BYTSCL( dd, min=min(dd[gd]), max=max(dd[gd]), top=252 )+1 ; 1 to 253
img[ bad ] = 254
erase

; rather than having a keyogram-ish type display with 12 (Jan -> Dec)
; rather than -90 -> 90, I'm making it 18 high. This way you can easily
; see an entire winter or summer season. Shift the addition by 1 year so that
; Dec of Year X moves up 1 pixel to Jan of Year X+1. Then, because
; we shifted, invalidate the shift wrap.
; img = [[img[*,0:*]],[SHIFT(img[*,0:5],-1)]]
; img[ n_elements(img[*,0])-1, 12:* ] = 254

!p.thick=3
!p.charthick=3
!p.charsize=1.0
!x.thick=3
!y.thick=3


position = [ 0.12, 0.55, 0.88, 0.925 ]

gd = where( s.annual ne 0 and s.annual ne 999.90 and s.years ne 0, ngd, complement=bad, ncomplement=nbad )
plotg, s.years[gd], s.annual[gd], $
       /yno, gval=224, $
       xrange=[1880,2010], /xst,psym=2, xminor=2, $
;       xticklen=-0.01, yticklen=-0.01, $
       ;xtitle='Time (Years)', $
       ytitle='Temperature (Degrees C)', $
       title=s.name, $
       position=position, thick=1, /nodata

;10 year smooth
ann = s.annual
if nbad ge 1 then ann[bad] = !values.f_nan
anns = smooth( ann, 10, /nan, /EDGE )
oplot, s.years[gd], anns[gd], thick=30, color=254

; overplot the asterisks
oplot, s.years[gd], s.annual[gd], psym=2, thick=1


; 1st order line fit to ALL data
gd = where( s.annual ne 999.90 AND s.years NE 0 )
line = indgen(n_elements(s.years))
r1 = poly_fit( line[gd], s.annual[gd], 1 )
fit1 = r1[0] + r1[1]*line
oplot, s.years[gd], fit1[gd],color=0

; 1st order line fit to last 50 years
gd = where( s.annual ne 999.90 and s.years ne 0 and s.years ge 2007-50, n )
if n ge 30 then begin
   line = indgen(n_elements(s.years))
   r1 = poly_fit( line[gd], s.annual[gd], 1 )
   fit1 = r1[0] + r1[1]*line
   oplot, s.years[gd], fit1[gd],color=150,thick=10
endif
; 1st order line fit to last 25 years
gd = where( s.annual ne 999.90 and s.years ne 0 and s.years ge 2007-25, n )
if n ge 15 then begin
   line = indgen(n_elements(s.years))
   r1 = poly_fit( line[gd], s.annual[gd], 1 )
   fit1 = r1[0] + r1[1]*line
   oplot, s.years[gd], fit1[gd],color=200,thick=10
endif
; 1st order line fit to last decade
gd = where( s.annual ne 999.90 and s.years ne 0 and s.years ge 2007-10, n )
if n ge 4 then begin
   line = indgen(n_elements(s.years))
   r1 = poly_fit( line[gd], s.annual[gd], 1 )
   fit1 = r1[0] + r1[1]*line
   oplot, s.years[gd], fit1[gd],color=250,thick=10
endif


position = [ position[0], 0.3, position[2], 0.475 ]
top = position[3]
imdisp, img, /noscale, aspect=0.2, /axis, /ystyle, order=1, $
        xrange=[1880,2010]-0.5, $
        xtickn=congrid([' '],60), $
        yticks=2, $
        ytickv=[0,5.5,11]+0.5, $
        ytickn=['Jan','Jun','Dec'], $
        yticklen=-0.01, xticklen=-0.1, $
;        ytitle='Time (Months)', $
        position=position, /usepos
;axis, /yaxis, yminor=1, yticklen=0.00001, yticks=1, ytickn=[' ',' ']

position[[1,3]] = position[[1,3]]-0.19
bot = position[1]
img = shift( img, 0, 6 )
img[ *, 0:5 ] = SHIFT( img[ *, 0:5 ], 1 )
img[ n_elements(img[*,0])-1, 0:5] = 254
imdisp, img, /noscale, aspect=0.2, /axis, /ystyle, order=1, $
        xrange=[1880,2010]-0.5, $
        xtickn=congrid([' '],60), $
        yticks=2, $
        ytickv=[0,5.5,11]+0.5, $
        ytickn=['July','Dec','Jun'], $
        yticklen=-0.01, xticklen=-0.1, $
;        xtitle='Time (Years)', $
;        ytitle='Time (Months)', $
        position=position, /usepos

bad = where( dd eq 0 or dd eq 999.90, nbad, complement=gd, ncomplement=ngd )
colorbar, bottom=1, $
          divisions=1, $
          maxrange=max(dd[gd]), $
          minrange=min(dd[gd]), $
          minor=1, $
          yticklen=1, $
          ncolors=252, $
          position=[ bot, 0.90, top, 0.92 ], $
          title='Temp (Deg. C)', $
          /right, /vertical
        


;gd = where( s.years ne 0 and s.annual ne 999.90 )
;axis, /yaxis, /save, yrange=minmax(s.annual[gd])+[-2,2], yminor=1, yticklen=-0.01
;oplot, s.years, s.annual, color=0, thick=2, psym=2

device, /close
;spawn, 'convert images/'+sname+'.ps ./images/'+sname+'.png'
set_plot,'x'
end






pro makeKML, giss, ghcn, data, meta
openw, lun, 'StationData_net.kml', /get_lun
printf, lun, '<?xml version="1.0" encoding="UTF-8"?>'
printf, lun, '<kml xmlns="http://earth.google.com/kml/2.1">'
printf, lun, '<Document>'

printf, lun, '<Folder>'
printf, lun, '<name>Station Data (2000 Anomaly)</name>'
for i=0,n_elements(giss)-1 do begin
;for i=0,100 do begin
    print, i

    giss_idx = i
    meta_idx =  where( giss[giss_idx].icc eq meta.icc AND $
                       giss[giss_idx].wmo eq meta.wmo AND $
                       giss[giss_idx].wmo_mod eq meta.wmo_mod, n )
    if n eq 0 then continue
    meta_idx = meta_idx[0]
    ghcn_idx =  where( giss[giss_idx].icc eq ghcn.icc AND $
                       giss[giss_idx].wmo eq ghcn.wmo AND $
                       giss[giss_idx].wmo_mod eq ghcn.wmo_mod, n )
    if n eq 0 then continue
    ghcn_idx = ghcn_idx[0]
    data_idx =  where( giss[giss_idx].icc eq data.icc AND $
                       giss[giss_idx].wmo eq data.wmo AND $
                       giss[giss_idx].wmo_mod eq data.wmo_mod, n )
    if n eq 0 then continue
    data_idx = data_idx[0]

;    if meta[meta_idx].nbaselineyr le 10 then continue
;    if meta[meta_idx].m25 eq 0 then continue

    if STRPOS( giss[giss_idx].name, "&" ) NE -1 then $
      giss[giss_idx].name = $
      STRJOIN(STRSPLIT(giss[giss_idx].name,'&',/EXTRACT),'and')
   
   printf, lun, '<Placemark>'
   printf, lun, '  <name>'+STRTRIM(giss[giss_idx].name,2)+'</name>'
   printf, lun, '  <Snippet></Snippet>'

   printf, lun, '  <description><![CDATA['
   description, lun, $
                giss[giss_idx].icc, $
                giss[giss_idx].wmo, $
                giss[giss_idx].wmo_mod, $
                giss[giss_idx].name
   printf, lun, '  ]]></description>'

   printf, lun, '  <Point>'
   printf, lun, '    <coordinates>'+ $
           STRTRIM(ghcn[ghcn_idx].lon,2)+','+STRTRIM(ghcn[ghcn_idx].lat,2)+ $
           '    </coordinates>'
   printf, lun, '  </Point>'


   icon, lun, giss[giss_idx], ghcn[ghcn_idx], meta[meta_idx], data[data_idx]

   printf, lun, '</Placemark>'
   printf, lun, ' '
endfor

printf, lun, '</Folder>'
printf, lun, '</Document>'
printf, lun, '</kml>'
free_lun, lun
end







pro description, lun,icc,wmo,wmo_mod,name
printf, lun, '    <center>'
printf, lun, '      <img src="http://edgcm.columbia.edu/~mankoff/StationData/images/'+ $
        icc+wmo+wmo_mod+'.ps.png">'
printf, lun, '    </center>'

spawn, 'grep '+icc+wmo+wmo_mod+' v2.temperature.inv.txt', out
printf, lun, '    <font size=-1>'
printf, lun, '      <pre>'
printf, lun, 'iccWMO_#... Name                              Lat     Lon Elev TEleP Pop Tp VLoCoAds------Vege------bi'
printf, lun, out
printf, lun, '<a href="http://edgcm.columbia.edu/~mankoff/StationData/v2.temperature.inv.txt">legend</a>'
printf, lun, '      </pre>'
printf, lun, '    </font>'

station = icc + wmo + wmo_mod
file = FINDFILE('data/'+station+'*.txt')
file = file[0]
file = strmid(file,strpos(file,'/',/reverse_search)+1,100) ; just the name
printf, lun, '    <font size=-1>' + $
        '<a href="http://edgcm.columbia.edu/~mankoff/StationData/data/'+ $
        file+'">station data file</a></font><br><br>'

printf, lun, '    <font size=-1>Produced by ' + $
        '<a href="mailto:mankoff@giss.nasa.gov">Ken Mankoff</a>'
printf, lun, 'for <a href="http://EdGCM.columbia.edu"</a>EdGCM</a>'
printf, lun, '    <br>'
printf, lun, '    Source: <a ' + $
        'href="http://data.giss.nasa.gov/gistemp/station_data/">NASA ' + $
        'GISS Temperature Station Data</a></font>'
end


pro icon, lun, giss, ghcn, meta, data

; load the colors. r,g,b are 16 long. 0 through 7 are 'cold' and
; 8 through 15 are 'warm'. White is not in the colorbar.
load_panoply,r,g,b

printf, lun, '<Style>'
printf, lun, '<IconStyle>'


; size is number of years of data
num_pts = n_elements(where( data.annual NE 999.90 AND data.years NE 0 ) )
size = KDM_RANGE( num_pts, 0.5, 1.5, $
                  datamin = 30, datamax = 2010-1880+30, $
                  underflow = 0.5, overflow = 1.5 )


; opacity is inverse proximity to city
opaque = KDM_RANGE( ghcn.pop, 128,255, $
                    datamin = 10, datamax = 20000, $
                    underflow = 128, overflow = 255 )
alpha = STRING( 255-opaque+128, FORMAT='(Z2.2)')



gd0 = where( data.years ge 1951 and $
             data.years le 1981 and $
             data.years ne 0 and $
             data.annual ne 999.9, n0 )
gd1 = where( data.years ge 2000 and $
             data.years ne 0 and $
             data.annual ne 999.9, n1 )
if n0 le 10 OR n1 lt 4 then begin
   bb = 'ff' & gg = 'ff' & rr = 'ff'
endif else begin
   base = mean( data.annual[gd0] )
   diff = mean( data.annual[gd1] )
   ctidx = KDM_RANGE( diff-base, 0,15, $
                      datamin=-3, datamax=3, $
                      underflow=0, overflow=15 )
   bb = b[ctidx] & gg = g[ctidx] & rr = r[ctidx]
   bb = string(bb,format='(Z2.2)')
   gg = string(gg,format='(Z2.2)')
   rr = string(rr,format='(Z2.2)')
;   if abs(diff-base) lt 0.1 then begin
;      bb = 'ff' & gg = 'ff' & rr = 'ff'
;   endif
endelse

printf, lun, '<scale>', + STRTRIM(size,2) + '</scale>'
printf, lun, '<color>'+alpha+bb+gg+rr+'</color>'

; Airport or not (use airplane icon)
printf, lun, '<Icon>'
;; Airplane icon shows up blue (Google Earth bug?), so only do pins...
;; printf, lun, ghcn.a eq 'A'?  $
;;         '<href>http://maps.google.com/mapfiles/kml/pal2/icon56.png</href>' : $
;;         '<href>http://maps.google.com/mapfiles/kml/paddle/wht-blank.png</href>'
printf, lun, '<href>http://maps.google.com/mapfiles/kml/paddle/wht-blank.png</href>'
printf, lun, '</Icon>'

printf, lun, '</IconStyle>'
printf, lun, '</Style>'
end

































pro write_vp
giss = load_giss_station_list(/restore)
ghcn = LOAD_GHCN_STATION_LIST(/restore)
data = load_station_data(/restore)
;meta = make_meta_data(/restore)

p_lut = ['R','S','U']
tp_lut = ['FL', 'HI', 'MT', 'MV']
v_lut = ['DE', 'FO', 'IC', 'MA', 'xx']
bi_lut = ['A','B','C']

openw, lun, 'vp.txt', /get_lun, width=512
printf, lun, '#lon lat ele tele pop co ds p tp v bi ' + $
        'nbaseyrs baseT recentT ' + $
        'trend_all, trend25'
for i=0, n_elements(giss)-1 do begin
   giss_idx = i
   ghcn_idx = where( giss[giss_idx].wmo eq ghcn.wmo AND $
                     giss[giss_idx].wmo_mod eq ghcn.wmo_mod AND $
                     giss[giss_idx].icc eq ghcn.icc, nghcn )
   if nghcn eq 0 then continue
   if nghcn ge 1 then ghcn_idx = ghcn_idx[0]
   data_idx = where( giss[giss_idx].wmo eq data.wmo AND $
                     giss[giss_idx].wmo_mod eq data.wmo_mod AND $
                     giss[giss_idx].icc eq data.icc, ndata )
   if ndata eq 0 then continue
   if ndata ge 1 then data_idx = data_idx[0]
;;     meta_idx = where( giss[giss_idx].wmo eq meta.wmo AND $
;;                       giss[giss_idx].wmo_mod eq meta.wmo_mod AND $
;;                       giss[giss_idx].icc eq meta.icc, nmeta )
;;     if nmeta eq 0 then continue
;;     if nmeta ge 1 then meta_idx = meta_idx[0]
   
; save GHCN everything and the data NumYears, baseline, and last decade
; ghcn.lat,lon,ele,tele,pop,co,ds

   p = where( ghcn[ghcn_idx].p eq p_lut )
   tp = where( ghcn[ghcn_idx].tp eq tp_lut )
   v = where( ghcn[ghcn_idx].v eq v_lut )
   bi = where( ghcn[ghcn_idx].bi eq bi_lut )
   
   ; record the number of baseline years
   nbaseyears = -999 & baselinetemp = -999
   baseyears = where( data[data_idx].years ge 1951 and $
                      data[data_idx].years le 1980 AND $
                      data[data_idx].years ne 0 AND $
                      data[data_idx].annual ne 999.9, nbaseyears )
   ; record the baseline temperature
   if nbaseyears gt 10 then $
      baselinetemp = mean( data[data_idx].annual[baseyears] )

   ; record the last 2k->present mean to compare with baseline
   recent_t = -999
   recent = where( data[data_idx].years ge 2000 and $
                   data[data_idx].years ne 0 and $
                   data[data_idx].annual ne 999.9, nrecent )
   if nrecent ge 4 then $
      recent_t = mean( data[data_idx].annual[recent] )
   
   ;; slope of 1st order line fit for all available years
   trend_all = -999
   gd = where( data[data_idx].years ne 0 AND $
               data[data_idx].annual ne 999.9, n )
   if n ge 2 then begin
      line = indgen(n_elements(data[data_idx].years))
      r1 = poly_fit( line[gd], data[data_idx].annual[gd], 1 )
      trend_all = r1[1]
   endif
   
   ;; trend for last 25 years
   trend25 = -999
   gd = where( data[data_idx].years ne 0 AND $
               data[data_idx].annual ne 999.9 AND $
               data[data_idx].years GE 1983, n )
   if n ge 2 then begin
      line = indgen(n_elements(data[data_idx].years))
      r1 = poly_fit( line[gd], data[data_idx].annual[gd], 1 )
      trend25 = r1[1]
   endif
   
   printf, lun, $
           ghcn[ghcn_idx].lon, $
           ghcn[ghcn_idx].lat, $
           ghcn[ghcn_idx].ele, $
           ghcn[ghcn_idx].tele, $
           ghcn[ghcn_idx].pop, $
           ghcn[ghcn_idx].co, $
           ghcn[ghcn_idx].ds, $
           p[0],tp[0],v[0],bi[0], $
           nbaseyears, baselinetemp, recent_t, $
           trend_all, trend25
   
endfor
free_lun, lun
end


pro stationdata

; get all the data
; 
; clean:
;; giss = LOAD_GISS_STATION_LIST()
;; ghcn = LOAD_GHCN_STATION_LIST()
;; data = load_station_data(giss,ghcn)
;; meta = make_meta_data(giss,ghcn,data)

; quick:
giss = LOAD_GISS_STATION_LIST(/restore)
ghcn = LOAD_GHCN_STATION_LIST(/restore)
data = load_station_data(/restore)
meta = make_meta_data(/restore)

;for i=0,100 do make_img, data[i]
for i=0,n_elements(data)-1 do make_img, data[i]
makeKML, giss, ghcn, data, meta

; $ for v in *.ps; do convert $v $v.png; done
; gzip the KML file
; upload KML_net, images, and if changed, KML and CGI

end







; color is temperature trend, and
; size is proximity to city (inverse)
; opacity is number of years of data

;; size is proximity to city (inverse)
;; size = KDM_RANGE( ghcn.pop, 0,1, $
;;                   datamin=10, datamax=20000, $
;;                   underflow=0, overflow=1) ; scale 0 to 1
;; size = 1.5-size ; invert. Cities are 0.5, rural is 1.5
;; printf, lun, '<scale>', + STRTRIM(size,2) + '</scale>'

;; ; opacity is number of years of data
;; num_pts = n_elements(where( data.annual NE 999.90 AND data.years NE 0 ) )
;; opaque = KDM_RANGE( num_pts, 128,255, datamin=25, datamax=2010-1880+30, $
;;                     underflow=128, overflow=255)
;; alpha = string(opaque,format='(Z2.2)')