; plotcrutsh.pro
; IDL program written by Tim Mitchell on 07.02.02
; last modified on 07.02.02

; ******************************************************************************
; unified entry

print, ""
print, "  > ***** PlotCRUtsH.pro: plots station locations *****"
print, ""

Black = 0 & Grey = 254 & White = 255
MissVal  = -999.0 & MissCol = White
ScaleSegN = 10 & BatchFileN = 0 & GivenYear = 0
InputText = "" & GivenFile = ""
RangeFract = 0.0D & ColSelect  = 0.0D

SeasonNames = strarr (4)
SeasonNames = ['winter (DJF)','spring (MAM)','summer (JJA)','autumn (SON)']

print, "  > Use standard (=0) or custom (=1) specifications ? "
read, QStanCust

; ******************************************************************************
; get basic stuff

print, "  > Enter the number of plots: "
read, PlotN

if            (PlotN EQ 1) then begin
  ColN = 1 & RowN = 1
endif else if (PlotN EQ 4) then begin
  ColN = 2 & RowN = 2
endif else if (PlotN EQ 6) then begin
  ColN = 3 & RowN = 2
endif else if (PlotN EQ 9) then begin
  ColN = 3 & RowN = 3
endif else begin
  print, "  > Specify the number of columns: "
  read, ColN  
  print, "  > Specify the number of rows: "
  read, RowN
endelse

;if (PlotN EQ 4) then begin
;  print, "  > Do the four plots correspond to the four seasons ? (0=no,1=yes)"
;  read, QSeasons
;endif else begin
;  QSeasons = 0
;endelse

QSeasons = 0

; ******************************************************************************
; get files and years

if (PlotN GT 1) then begin
  print, "  > Do we vary by .cts file (=0) or by year (=1) ?
  read, QFileYear
endif else begin
  QFileYear = -1
endelse

FilePaths = strarr (PlotN)
if (QFileYear EQ 0) then begin
  for XPlot = 0, (PlotN-1) do begin
    print, "  > Enter the .cts or .hdr file to load for plot: ", (XPlot+1)
    read, GivenFile
    FilePaths[XPlot] = GivenFile
  endfor
endif else begin
  print, "  > Enter the .cts or .hdr file to load: "
  read, GivenFile
  for XPlot = 0, (PlotN-1) do begin
    FilePaths[XPlot] = GivenFile
  endfor
endelse

FileYears = intarr (PlotN)
if (QFileYear EQ 1) then begin
  print, "  > Enter a sequence (start,end) or individual yrs (-999,-999): "
  read, SeqBeg,SeqEnd
  
  if (SeqBeg Eq -999 OR SeqEnd EQ -999) then begin
   for XPlot = 0, (PlotN-1) do begin
    print, "  > Enter the year AD (-999=any) for plot: ", (XPlot+1)
    read, GivenYear
    FileYears[XPlot] = GivenYear
   endfor
  endif else begin
   for XPlot = 0, (PlotN-1) do begin
    FileYears[XPlot] = SeqBeg + XPlot
   endfor
  endelse
endif else begin
  print, "  > Enter the year AD (-999=any): "
  read, GivenYear
  for XPlot = 0, (PlotN-1) do begin
    FileYears[XPlot] = GivenYear
  endfor
endelse

GetViewBounds, ViewN, ViewName, ViewKeyWords, ViewBounds, ViewFullName=ViewFullName

LoadAnyCT, 50

; GetScales, ScaleN, ScaleName, ScaleColTab, ScaleSeg, ScaleLimits 

; ******************************************************************************
; plot choices

print, "  > Do we label the plots ? (0=no,1=yes)" 
read, QPlotName
 
if (QPlotName EQ 1 AND QSeasons EQ 0) then begin
   if (QFileYear EQ 1) then begin
     print,"  > Label with yrs (=1), yrs+tots (=2) ? (manual=0)" 
     read,QLabelYear
   endif else begin
     QLabelYear = 0
   endelse
   
   PlotName = strarr (PlotN)
   if (QLabelYear EQ 0) then begin
     print, "  > Enter the title for each plot in turn: "
     for XPlot = 0, (PlotN-1) do begin
      read, InputText
      PlotName[XPlot] = InputText
      GiveSymbol, PlotName[XPlot]
     endfor
   endif else begin
     for XPlot = 0, (PlotN-1) do begin
      PlotName[XPlot] = string(FileYears[XPlot])
     endfor
   endelse
endif

;if (PlotN NE 1) then begin
;  print, "  > Do the plots have the same scale ? (1=yes,2=no)"
;  read, QSameScale
;endif else begin
;  QSameScale = 1
;endelse
;if (QSameScale EQ 1) then begin
;  Info = strarr (10) & Info [*] = "missing"
;  GrabInfo, "@@@@@@@@", "@@@@@@@@", Info, -999, 0
;  
;  SeasName = "missing"
;  SearchSeas, FilePaths(0), FileInfos(0), SeasName, SeasIndex=SeasIndex
;endif
ContourFactor = 5
Info = strarr (10) & Info [*] = "missing"
Info[0]="elevation" & Info[1]="mean" & Info[3]="metres"
LoadVariCT, Info, ContourFactor, DataBounds, ContourN, TickVals, TickFormat
ColorBounds = [1,253]
GetContours, DataBounds, ColorBounds, Contours, ContourColors, ContourN
QSameScale = 1

QView = MissVal & InputInt = 0
while (QView LT 0) do begin
     print, "  > Identify the region to plot (-1=list): "
     read, InputInt
    
     if (InputInt EQ -1) then begin
      for XView = 0, (ViewN-1) do begin
        print, "  > ", XView, " : ", ViewName[XView], format='(a4,i4,a3,a10)'
      endfor
     endif else begin
      if (InputInt GE 0 AND InputInt LT ViewN) then QView = InputInt
     endelse
endwhile
Limits = [ViewBounds[QView,2],ViewBounds[QView,0],ViewBounds[QView,3],ViewBounds[QView,1]]

print, "  > Label with station names ? (0=no,1=yes)"
read, QLabelName
       
;print, "  > Plot onto map (=0) or grid (=1) ? "
;read, QMapGrid
QMapGrid = 0

;print, "  > Plot contours (=0) or individual grid cells (=1) ?"
;read, QContourCell

;print, "  > Select the contour density (1=min,5=max): "
;read, ContourFactor

if (QMapGrid EQ 0) then begin
;  Indent = 0.02
  print, "  > Plot grid ? (0=no,1=yes)"
  read, QGrid
;  if (QGrid EQ 1) then Indent = 0.05
  
  print, "  > Plot outlines ? (-2=IDLcoasts,-1=select,0=no,>0=predefined)"
  read, QOutline
  if (QOutLine EQ -1 OR QOutLine GT 0) then SelectPoly, QOutline
endif else begin
  QGrid = -999 & QOutline = -999
;  Indent = 0.02
endelse
Indent = 0.02

print, "  > Side bar: 0=none,1=all,2=attrib,3=info (scale is plotted)"
read, QBarInfo

print, "  > Plot shape: standard (=0) or fit to the region plotted (=1) ?"
read, QStanFit

;if (QStanFit EQ 1 AND QOneGrid EQ 2) then begin
;  print, "  > Plot shape: enter LongN,LatN: "
;  read, LongN,LatN
;endif
LongN = Limits[3]-Limits[1] & LatN = Limits[2]-Limits[0]

print, "  > Plot  size factor (standard=1.0)  ?"
read, PlotSizeFactor
PlotSizeFactor = PlotSizeFactor / 200.0

LabelSizeFactor = 1.0
print, "  > Label size factor (standard=1.0)  ?"
read, LabelSizeFactor

TextSize = 0.8 & TextThick = 2.0

;print, "  > Plot to screen (=21) or .eps (=22) ?"
;read, Choice
Choice = 22

; ******************************************************************************
; determine plot size and identify plot segments
; note that when multiple grids are used in a single execution, this will go wonky

if (QBarInfo EQ 0) then begin 				; this allows width for just the scale
  BarPixels = 2.0
endif else begin					; this allows width for the logo and info as well
  BarPixels = 4.0
endelse

if (QPlotName EQ 0) then begin 				; this allows height for the label, if required
  NamePixels = 0.0
endif else begin					
  NamePixels = 1.0
endelse

DataWidth = Limits[3] - Limits[1]
DataHeight= Limits[2] - Limits[0]

ShrinkRatio = sqrt((DataWidth*DataHeight)/150.0)
DataWidth = 200.0 * DataWidth / ShrinkRatio
DataHeight= 200.0 * DataHeight/ ShrinkRatio

if (QStanFit EQ 0) then begin				; this makes any PS plot have a size 16cm * 12cm
							;   so fits into Word doc without any resizing
  XPixels = (16.0*PlotSizeFactor) & YPixels = (12.0*PlotSizeFactor)  	
endif else begin					; this sizes the plot according to the .glo sizes
					  		;   so minimises white space in the plot
  YPixels = float(RowN) * PlotSizeFactor * (DataHeight + NamePixels)
  
  if (QSameScale EQ 2) then begin
    XPixels = PlotSizeFactor * ((float(ColN)*DataWidth) + (1.0*BarPixels))
  endif else begin
    XPixels = PlotSizeFactor * ((float(ColN)*DataWidth) + (float(ColN)*BarPixels))
  endelse
endelse

BarWidth   = BarPixels / XPixels			; this converts the key width into normalised units
NameHeight = LabelSizeFactor * NamePixels / YPixels	; ...and the name height likewise

							; this identifies the normalised key boundaries
if (QSameScale EQ 2) then begin				; here, one scale per .glo
  GetCorners, [0.0,0.0,1.0,1.0], PlotCorners, ColN, RowN, 0, Indent=Indent
  ScaleCorners = PlotCorners
  PlotCorners  [*,2] = PlotCorners  [*,2] - BarWidth
  ScaleCorners [*,0] = ScaleCorners [*,2] - BarWidth
endif else begin					; here one scale for the entire plot
  GetCorners, [0.0,0.0,(1.0-BarWidth),1.0], PlotCorners,  ColN, RowN, 0, Indent=Indent
  GetCorners, [(1.0-BarWidth),0.0,1.0,1.0], ScaleCorners, 1,       1, 0, Indent=Indent
endelse

if (QPlotName EQ 1) then begin				; this provides space for a label for each .glo
  LabelCorners = PlotCorners
  
  PlotCorners  [*,3] = PlotCorners  [*,3] - NameHeight
  LabelCorners [*,1] = LabelCorners [*,3] - NameHeight

  if (QSameScale EQ 2) then ScaleCorners [*,3] = ScaleCorners [*,3] - NameHeight
endif else begin
  LabelCorners = PlotCorners * 0.0
endelse

for XPlot = 0, (PlotN-1) do begin			; this ensures each .glo has correct width/length
  XGloPixels = (PlotCorners[XPlot,2] - PlotCorners[XPlot,0]) * XPixels
  YGloPixels = (PlotCorners[XPlot,3] - PlotCorners[XPlot,1]) * YPixels
  
  if (YGloPixels/XGloPixels LT float(LatN)/float(LongN)) then begin
    XNewPixels = YGloPixels * (float(LongN)/float(LatN))
    PlotCorners[XPlot,0] = PlotCorners[XPlot,0] + ((XGloPixels-XNewPixels)/(2.0*XPixels))
    PlotCorners[XPlot,2] = PlotCorners[XPlot,2] - ((XGloPixels-XNewPixels)/(2.0*XPixels))
  endif else begin
    YNewPixels = XGloPixels * (float(LatN)/float(LongN))
    PlotCorners[XPlot,1] = PlotCorners[XPlot,1] + ((YGloPixels-YNewPixels)/(2.0*YPixels))
    PlotCorners[XPlot,3] = PlotCorners[XPlot,3] - ((YGloPixels-YNewPixels)/(2.0*YPixels))
  endelse
endfor

SymVec = findgen(16) * (!pi*2/16.)
UserSym, cos(SymVec), sin(SymVec), /fill

; ******************************************************************************
; plot by plot

SeasName = "" & SeasClash = 0

if (Choice EQ 21) then begin
  set_plot, 'X'	
  if (!D.Window EQ -1) then window, /free, xsize=(XPixels*32.0), ysize=(YPixels*32.0), $
  		xpos=(1000-(XPixels*32.0)), title="box plot"
endif else begin
  SavePath = ""
  print, "  > Enter the .eps or .cgm file to save: "
  read, SavePath

  EPSbeg = strpos(SavePath,".eps")
  CGMbeg = strpos(SavePath,".cgm")
  
  if (EPSbeg ne -1) then begin
    Set_Plot, 'ps', /copy
    device, filename=SavePath, bits_per_pixel=8, xsize=XPixels, ysize=YPixels, /Color, /encapsulated
  endif else if (CGMbeg ne -1) then begin
    Set_Plot, 'cgm', /copy
    device, filename=SavePath
  endif
endelse

for XPlot = 0, (PlotN-1) do begin
  if (QSeasons NE 1) then begin			; these are not the four seasons
    print, "  > Plot: ", (XPlot+1)
  endif else begin				; these are the four seasons
    print, "  > Plot: ", SeasonNames (XPlot)
  endelse
  
  ExeRef = ((LabelCorners[XPlot,2] - LabelCorners[XPlot,0]) / 2.0) + LabelCorners[XPlot,0]
  WyeRef = ((LabelCorners[XPlot,3] - LabelCorners[XPlot,1]) / 3.0) + LabelCorners[XPlot,1]
  
  if (XPlot EQ 0) then begin
      LoadCRUtsHeads, Code,LatFlt,LonFlt,ElvFlt,YearAD0,YearAD1,NameStn,NameCty, CallFile=FilePaths(XPlot)
  endif else begin
      if (FilePaths(XPlot) NE FilePaths(XPlot-1)) then $
      		LoadCRUtsHeads, Code,LatFlt,LonFlt,ElvFlt,YearAD0,YearAD1,NameStn,NameCty, CallFile=FilePaths(XPlot)
  endelse

; ##### do plot #####
  
  if (QMapGrid EQ 0) then begin
    if (ViewKeyWords[QView,0] EQ 'Lamb') then Map_set, 0, 0, /lambert, color=Black, /noerase, $
  	position=PlotCorners[XPlot,*], /noborder, limit=Limits 
    if (ViewKeyWords[QView,0] EQ 'Cyli') then Map_set, 0, 0, /cylindrical, color=Black, /noerase, $
  	position=PlotCorners[XPlot,*], /noborder, limit=Limits
    if (ViewKeyWords[QView,0] EQ 'Hamm') then Map_set, 0, 0, /hammer, color=Black, /noerase, $
  	position=PlotCorners[XPlot,*], /noborder, limit=Limits
  endif
  
  if (FileYears[XPlot] NE -999) then begin
    Combo = (FileYears[XPlot]-YearAD1(*)) * (YearAD0(*)-FileYears[XPlot])
  endif else begin
    Combo = YearAD0(*)
  endelse

  Legit = where(Combo NE -999, LegitN)
  for XLegit = 0, (LegitN-1) do begin
      if (Combo[XLegit] GE 0) then begin
        if (LonFlt[XLegit] LT Limits[1] OR LonFlt[XLegit] GT Limits[3] OR LatFlt[XLegit] LT Limits[0] $
        		OR LatFlt[XLegit] GT Limits[2]) then Combo[XLegit] = -999
      endif
  endfor
    
  Include = where(Combo GE 0, IncludeN)

  LonInclude = LonFlt(Include)
  LatInclude = LatFlt(Include)
  ElvInclude = ElvFlt(Include)
  NamInclude = NameStn(Include)
  
  for XInclude = 0, (IncludeN-1) do begin
    if (ElvInclude(XInclude) NE MissVal) then begin
      if (ElvInclude(XInclude) LT DataBounds[0]) then begin
        ElvInclude(XInclude) = ColorBounds[0]
      endif else begin
        if (ElvInclude(XInclude) GT DataBounds[1]) then begin
        	ElvInclude(XInclude) = ColorBounds[1]
        endif else begin
        	RangeFract = (ElvInclude(XInclude)-DataBounds[0])/(DataBounds[1]-DataBounds[0])
        	ColSelect  = (RangeFract*float(DataBounds[1]-DataBounds[0])) + float(DataBounds[0])
        	ElvInclude(XInclude) = round (ColSelect)
        endelse
      endelse 
    endif else begin
      ElvInclude(XInclude) = Black
    endelse
  endfor  
    
  print, "  > Total stations plotted: ", IncludeN
  if (QLabelYear EQ 2) then PlotName[XPlot]=strtrim(PlotName[XPlot],2)+': '+string(IncludeN)+' stations'
  if (QPlotName EQ 1) then xyouts, ExeRef,WyeRef,PlotName[XPlot],font=5,charsize=(1.2*LabelSizeFactor),charthick=2.0,/normal,alignment=0.5
  
  LonOffset = (Limits[3]-Limits[1])/100.0
  LatOffset = (Limits[2]-Limits[0])/100.0
  plots,LonInclude,LatInclude,PSym=8,SymSize=0.3,Color=ElvInclude
  if (QLabelName EQ 1) then xyouts,(LonInclude+LonOffset),LatInclude,NamInclude,Color=Black,CharSize=0.5
  
  if (QOutline NE 0 AND QMapGrid EQ 0) then begin
   if (QOutline EQ -2) then begin
    if (ViewKeyWords[QView,0] EQ 'Lamb') then Map_set, 0, 0, /lambert, color=Black, /noerase, $
  	position=PlotCorners[XPlot,*], /continents, /noborder, /hires, limit=Limits 
    if (ViewKeyWords[QView,0] EQ 'Cyli') then Map_set, 0, 0, /cylindrical, color=Black, /noerase, $
  	position=PlotCorners[XPlot,*], /continents, /noborder, /hires, limit=Limits
    if (ViewKeyWords[QView,0] EQ 'Hamm') then Map_set, 0, 0, /hammer, color=Black, /noerase, $
   	position=PlotCorners[XPlot,*], /continents, /noborder, /hires, limit=Limits
   endif else begin
    PlotPoly, QOutline, linestyle=0, color=Black
   endelse
  endif
  
  if (QGrid EQ 1 AND QMapGrid EQ 0) then begin
    if (ViewKeyWords[QView,0] EQ 'Lamb') then Map_set, 0, 0, /lambert, color=Black, /noerase, $
  	position=PlotCorners[XPlot,*], /grid, /label, glinestyle=0, /noborder, limit=Limits, $
  	LatLab=(Limits[1]+LonOffset), LonLab=(Limits[0]+LatOffset), latalign=0,lonalign=0,charsize=0.5
    if (ViewKeyWords[QView,0] EQ 'Cyli') then Map_set, 0, 0, /cylindrical, color=Black, /noerase, $
  	position=PlotCorners[XPlot,*], /grid, /label, glinestyle=0, /noborder, limit=Limits, $
  	LatLab=(Limits[1]+LonOffset), LonLab=(Limits[0]+LatOffset), latalign=0,lonalign=0,charsize=0.5
    if (ViewKeyWords[QView,0] EQ 'Hamm') then Map_set, 0, 0, /hammer, color=Black, /noerase, $
   	position=PlotCorners[XPlot,*], /grid, /label, glinestyle=0, /noborder, limit=Limits, $
  	LatLab=(Limits[1]+LonOffset), LonLab=(Limits[0]+LatOffset), latalign=0,lonalign=0,charsize=0.5
  endif
  
  if (QSameScale EQ 2) then DrawConScale,ScaleCorners[XPlot,*],Contours,ContourColors,$
    	TickVals,TickFormat,TextSize,TextThick,Info,XPixels,YPixels,QBarInfo
endfor

if (QSameScale EQ 1) then DrawConScale,ScaleCorners[0,*],Contours,ContourColors,$
    	TickVals,TickFormat,TextSize,TextThick,Info,XPixels,YPixels,QBarInfo

; ******************************************************************************
; end plot

if (Choice ne 21) then device, /close

set_plot, 'X'	

; ******************************************************************************
; end program

print, ""

end