;
; Makes plots of tree site locations for my webpage.
; Use ps2gif <fn> 0.5 for big figures
; Use ps2gif <fn> 0.1 for thumbnails
;
;iver=0   ; 0=use regional lists, 1=use overall list
;
; Get sites used in Briffa et al. (2002) Holocene paper, to see which ones are
; the same.  Of course they should all be the same, but since there I dealt
; with standardised chronologies provided by Schweingruber, while here I have
; the raw core measurements, there's no guarantee they were all used in both
; data sets!
;
; Get site locations
;
restore,filename='/cru/u2/f055/tree6/allmxd.idlsave'
;  nchron,idno,idname,location,country,tree,yrstart,yrend,statlat,statlon,$
;  mxd,fraction,timey,nyr
;
hnchron=nchron
hidname=idname
hlocation=string( ( byte(location) < 126B ) > 32B )
htree=strcompress(tree,/remove_all)
hyrst=yrstart
hyren=yrend
hlat=statlat
hlon=statlon
;
; Get site selection list
;
nreg=9
openr,1,'ninelist.scr'
headst=''
readf,1,headst
for ireg = 0 , nreg-1 do begin
 readf,1,nloc,format='(I4)'
 loclist=fltarr(2,nloc)
 readf,1,loclist,format='(I5,1X,I6)'
 if ireg eq 0 then begin
  locx=reform(loclist[1,*]/100.)
  locy=reform(loclist[0,*]/100.)
  nsel=nloc
 endif else begin
  locx=[locx,reform(loclist[1,*]/100.)]
  locy=[locy,reform(loclist[0,*]/100.)]
  nsel=nsel+nloc
 endelse
endfor
close,1
;
; Get site locations
;
; First read header lines of *all* core data
;
fndat='mxd.allsites.all.dat'
print,fndat
;
ncore=0
idno1=''
meas1=''
proc1=''
tree1=''
idname1=''
location1=''
;
openr,1,fndat
;
while not eof(1) do begin
;
headst=''
readf,1,idno1,yrstart1,yrend1,meas1,proc1,tree1,pith1,statlat1,$
 statlon1,elev1,idname1,location1,mpr1,$
 format='(A6,2I4,2X,3A4,I4,I5,I6,I4,A8,A25,I4)'
yrstart1=fix(yrstart1)
yrend1=fix(yrend1)
nline=(yrend1/10)-(yrstart1/10)+1
bodyst=strarr(nline)
readf,1,bodyst
;
; Only add these site details to the list if it's different to the
; immediately preceding core
;
addit=0
if ncore eq 0 then begin
 addit=1
endif else begin
 if (idname1 ne idname[ncore-1]) or (tree1 ne tree[ncore-1]) or $
  (statlat1 ne statlat[ncore-1]) or (statlon1 ne statlon[ncore-1]) then begin
  addit=1
 endif else begin
  ; Although it is from the same site as the last core and so we don't need to
  ; add it to the list, we might want to update the start and end years of the
  ; overall site chronology, in case the new core extends earlier or later than
  ; the other cores
  yrstart[ncore-1]=yrstart[ncore-1] < yrstart1
  yrend[ncore-1]=yrend[ncore-1] > yrend1
 endelse
endelse
;
if addit eq 1 then begin
 if ncore eq 0 then begin
  idno=[idno1]
  yrstart=[yrstart1]
  yrend=[yrend1]
  tree=[tree1]
  statlat=[statlat1]
  statlon=[statlon1]
  idname=[idname1]
  location=[location1]
 endif else begin
  idno=[idno,idno1]
  yrstart=[yrstart,yrstart1]
  yrend=[yrend,yrend1]
  tree=[tree,tree1]
  statlat=[statlat,statlat1]
  statlon=[statlon,statlon1]
  idname=[idname,idname1]
  location=[location,location1]
 endelse
 ncore=ncore+1
endif
;
endwhile
;
close,1
;
print,'NUMBER OF SITES:',ncore
ploton,1
;
; Get rid of control codes from location names!
;
location=string( ( byte(location) < 126B ) > 32B )
;
statlat=float(statlat)/100.
statlon=float(statlon)/100.
;
masklist=fltarr(ncore)*!values.f_nan
for icore = 0 , ncore-1 do begin
 gotit=where( (locx eq statlon[icore]) and (locy eq statlat[icore]),ngot)
 case ngot of
  0: masklist[icore]=0
  else: masklist[icore]=1
 endcase
endfor
print,'NUMBER AFTER SELECTION:',total(masklist)
;
loadct,39
multi_plot,nrow=2,ncol=1,layout='large'
if !d.name eq 'X' then begin
  window,ysize=850
  fac=1.
endif else begin
  fac=2.
endelse
def_1color,10,color='vvlgreen'
;
map=def_map(/npolar)  &  map.limit(0)=25.
labels=def_labels(/off)  &  labels.gridon=0  &  labels.dlon=30.
labels.dlat=10.
coast=def_coast(/get_device)
coast.double=0
coast.fill=1
coast.fillcolor=10
coast.thick=[0.5,0.5]
;
inter_boxfd,/nodata,map=map,$
  labels=labels,coast=coast
;
!p.ticklen=0.03
lon_polar,map=map,[-180,-150.,-120.,-90.,-60.,-30.,30.,60.,120.,150.],/dotick
;
; Overlay regions
;
map_plots,[-15.,-15.],[90.,35.],thick=3
map_plots,[61.,61.],[90.,60.],thick=3
map_plots,[128.,128.],[90.,60.],thick=3
map_plots,[-170.,-170.],[90.,50.],thick=3
map_plots,[-125.,-125.],[90.,60.],thick=3
map_plots,[-45.,-45.],[90.,40.],thick=3
map_plots,[30.,30.],[60.,35.],thick=3
map_plots,[115.,115.],[60.,25.],thick=3
map_plots,[60.,60.],[40.,25.],thick=3
map_plots,[-135.,-135.],[50.,30.],thick=3
map_plots,[-100.,-100.],[50.,30.],thick=3
map_plots,[-105.,-105.],[55.,50.],thick=3
map_plots,[-120.,-120.],[60.,55.],thick=3
;
y1=[60.,50.,30.,60.,55.,40.,53.,35.,60.,40.,25.]
x1=[-180.,-170.,-135.,-125.,-120.,-100.,-15.,-15.,30.,30.,60.]
x2=[-170.,-100.,-100.,-120.,-105.,-45.,30.,30.,180.,115.,115.]
for i = 0 , n_elements(y1)-1 do begin
  n=x2(i)-x1(i)+1
  xval=findgen(n)+x1(i)
  yval=replicate(y1(i),n)
  map_plots,xval,yval,thick=3
endfor
;
regname=['NEUR','SEUR','NSIB','ESIB','CAS','TIBP',$
  'WNA','NWNA','ECCA']
reglat=[  73,  35,  85,  81,  53,  37,  30,  67,  62]
reglon=[ -10, -19,  86,-176,  38,  70,-135,-129, -72]
for i = 0 , n_elements(regname)-1 do $
  xyouts,reglon(i),reglat(i),regname(i),charsize=0.7
;
xcir=findgen(361)-180.
ycir=replicate(map.limit[0],361)
map_plots,xcir,ycir
;
; Output table of sites used
;
 ns=ncore
nst=txt(ns)
;
openw,1,'b01abd_site.txt'
printf,1,'Table of tree-ring sites used in'
printf,1
printf,1,'Briffa KR et al. (2001) Low-frequency temperature variations from a northern tree-ring-density network. J. Geophysical Research 106, 2929-2941.'
printf,1
printf,1,'Number of sites: '+nst
printf,1
printf,1,'                      Name     Abbr Genus Longitude  Latitude Start   End'
for is = 0 , ns-1 do printf,1,strtrim(location[is],2),$
 strtrim(idname[is],2),strtrim(tree[is],2),statlon[is],statlat[is],$
 yrstart[is],yrend[is],$
 format='(A26,A9,A6,2F10.2,2I6)'
close,1
;
x=statlon
y=statlat
cpl_usersym,/circle,/fill
plots,x,y,psym=8,color=!p.background,symsize=1.05,noclip=0
cpl_usersym,/circle
plots,x,y,psym=8,thick=1,symsize=1.05,noclip=0
print,'Number of sites:',n_elements(x)
;
plotoff
;
; Now go through list from the Holocene paper and cross check for differences
; NB hugershoff chronologies had lat/lon in degrees&minutes, while core data
; have it in degrees&hundredths, so location will only match to 2 d.p.
;
openw,1,'comparesites_abd_holocene.txt'
printf,1,'Those in B01abd that are not in B02sens'
mslist=fltarr(hnchron)*!values.f_nan
nkeep=0
for is = 0 , ns-1 do begin
 ;
 matchlist=where( (abs(statlon[is]-hlon) lt 0.01) and (abs(statlat[is]-hlat) lt 0.01) and $
  (yrstart[is] eq hyrst) and (yrend[is] eq hyren) and (tree[is] eq htree) , nmatch)
 ;
 if nmatch eq 0 then begin
  printf,1,strtrim(location[is],2),$
   strtrim(idname[is],2),strtrim(tree[is],2),statlon[is],statlat[is],$
   yrstart[is],yrend[is],$
   format='(A26,A9,A6,2F10.2,2I6)'
  nkeep=nkeep+1
 endif else begin
  mslist[matchlist]=1
  if nmatch gt 1 then begin
    printf,1,'MULTIPLE MATCHES!'
   for js = 0 , nmatch-1 do begin
    ks=matchlist[js]
  printf,1,strtrim(hlocation[ks],2),$
   strtrim(hidname[ks],2),strtrim(htree[ks],2),hlon[ks],hlat[ks],$
   hyrst[ks],hyren[ks],$
   format='(A26,A9,A6,2F10.2,2I6)'
   endfor
    printf,1,'MULTIPLE MATCHES!'
   endif
 endelse
 ;
endfor
printf,1,'Total=',nkeep
printf,1
printf,1,'Those in B02sens that are not in B01abd'
kl=where(finite(mslist) eq 0,nkeep)
for ij = 0 , nkeep-1 do begin
  is=kl[ij]
  printf,1,strtrim(hlocation[is],2),$
   strtrim(hidname[is],2),strtrim(htree[is],2),hlon[is],hlat[is],$
   hyrst[is],hyren[is],$
   format='(A26,A9,A6,2F10.2,2I6)'
endfor
printf,1,'Total=',nkeep
close,1
;
spawn,'cp idl.ps1 b01abd_map.ps'
spawn,'mv idl.ps1 b01abd_map_thumb.ps'
spawn,'ps2gif b01abd_map 0.5'
spawn,'ps2gif b01abd_map_thumb 0.1'
;
end