;+
; NAME:
;       PrecipMap
;
; PURPOSE:
;
;       This is a program that demonstrates how to place an IDL map projection
;       onto an image that is already in a map projection space. It uses a NOAA
;       precipitation image that is in a polar stereographic map projection, and
;       for which the latitudes and longitudes of its four corners are known.
;
;       For additional details, see http://www.idlcoyote.com/map_tips/precipmap.html.
;
; AUTHOR:
;
;       FANNING SOFTWARE CONSULTING
;       David Fanning, Ph.D.
;       1645 Sheely Drive
;       Fort Collins, CO 80526 USA
;       Phone: 970-221-0438
;       E-mail: david@idlcoyote.com
;       Coyote's Guide to IDL Programming: http://www.idlcoyote.com/
;
; CATEGORY:
;
;       Map Projection
;
; CALLING SEQUENCE:
;
;       IDL> PrecipMap, filename
;
; INPUTS:
;
;      filename:   The name of the precipitation file. For demo, download
;                  ST4.2005010112.24h.bin from http://www.idlcoyote.com/data/ST4.2005010112.24h.bin.
;                  
; KEYWORDS:
; 
;      DATA:   Set this keyword to a named variable that on output will contain the scaled data.
;      
;      PALETTE:    Set this keyword to a named variable that on output will contain the color 
;                  palette used to display the data.
;                   
;
; RESTRICTIONS:
;
;     Requires files from the Coyote Library:
;
;     http://www.idlcoyote.com/documents/programs.html
;
; MODIFICATION HISTORY:
;
;  Written by David W. Fanning, 28 April 2006 from code and discussion supplied
;       by James Kuyper in the IDL newsgroup.
;  Renamed Colorbar procedure to cgColorbar to avoid conflict with IDL 8 Colorbar function.
;        26 September 2010. DWF.
;  Got the polar stereo map projection correct. 5 September 2011. DWF.
;  Added DATA, MAP, and PALETTE output keywords and updated to use more modern Coyote
;        Library mapping routines. 2 November 2012. DWF.
;-
;
;******************************************************************************************;
;  Copyright (c) 2008, by Fanning Software Consulting, Inc.                                ;
;  All rights reserved.                                                                    ;
;                                                                                          ;
;  Redistribution and use in source and binary forms, with or without                      ;
;  modification, are permitted provided that the following conditions are met:             ;
;                                                                                          ;
;      * Redistributions of source code must retain the above copyright                    ;
;        notice, this list of conditions and the following disclaimer.                     ;
;      * Redistributions in binary form must reproduce the above copyright                 ;
;        notice, this list of conditions and the following disclaimer in the               ;
;        documentation and/or other materials provided with the distribution.              ;
;      * Neither the name of Fanning Software Consulting, Inc. nor the names of its        ;
;        contributors may be used to endorse or promote products derived from this         ;
;        software without specific prior written permission.                               ;
;                                                                                          ;
;  THIS SOFTWARE IS PROVIDED BY FANNING SOFTWARE CONSULTING, INC. ''AS IS'' AND ANY        ;
;  EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES    ;
;  OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT     ;
;  SHALL FANNING SOFTWARE CONSULTING, INC. BE LIABLE FOR ANY DIRECT, INDIRECT,             ;
;  INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED    ;
;  TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;         ;
;  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND             ;
;  ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT              ;
;  (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS           ;
;  SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.                            ;
;******************************************************************************************;

FUNCTION PrecipMap_Annotate, axis, index, value

   ; This is a function for annotating the colorbar.

   plevels = Indgen(15)
   levels = [0, 2, 5, 10, 15, 20, 25, 50, 75, 100, 125, 150, 200, 250]
   labels = StrTrim(levels, 2)
   text = [labels[0:12], labels[13] + '+', '']
   selection = Where(plevels EQ value)
   RETURN, (text[selection])[0]

END ;----------------------------------------------------------------------



PRO PrecipMap, filename, DATA=data, MAP=map, Palette=palette

   ; Error handling.
   Catch, theError
   IF theError NE 0 THEN BEGIN
      Catch, /Cancel
      void = cgErrorMsg()
      RETURN
   ENDIF

   ; Need a filename?
   IF N_Elements(filename) EQ 0 THEN BEGIN
      filename = Dialog_Pickfile(Filter='*.bin') ; Demo file.
      IF filename EQ "" THEN RETURN
   ENDIF
   OPENR, lun, filename, /Get_Lun
   image = FltArr(1121, 881)
   ReadU, lun, image
   Free_Lun, lun

   ; Set up 15 colors, roughly equvilent to NOAA precipitation colors.
   ; Use black for missing values.
   colors = ['dark green', 'lime green', 'light sea green', 'yellow', 'khaki', $
             'dark goldenrod', 'light salmon', 'orange', 'red', 'sienna', 'purple', $
             'orchid', 'thistle', 'sky blue', 'black']
   TVLCT, cgColor(colors, /Triple), 1

   ; Process image data. Missing values = 9.999e20.
   s = Size(image, /Dimensions)
   scaledImage = BytArr(s[0],s[1])
   levels = [0, 2, 5, 10, 15, 20, 25, 50, 75, 100, 125, 150, 200, 250]
   FOR j=1,13 DO BEGIN
      index = Where(image GE levels[j-1] AND image LT levels[j], count)
      IF count GT 0 THEN scaledImage[index] = Byte(j)
   ENDFOR
   index = Where(image GT 250 AND image LT 9.999e20, count)
   IF count GT 0 THEN scaledImage[index] = 14B
   missing = Where(image EQ 9.999e20, missing_count)
   IF missing_count GT 0 THEN scaledImage[missing] = 15B ; Will be black.

   ; Set up the map structure. This image is a stereographic image.
   mapCoord = Obj_New('cgMap', "Polar Stereographic", CENTER_LONGITUDE=-105, CENTER_LATITUDE=70)

   ; Set up latitutes and longitudes at the corners of the image. (ll, ul, ur, lr)
   longitude = [-119.023D, -134.039, -59.959, -80.7500]
   latitude =  [  23.117D,   53.509,  45.619,  19.8057]

   ; Project those lat/lon points into UV space.
   uv = mapCoord -> Forward(longitude, latitude)
   
   ; To set up map projection space, we need values at left, top, right, and bottom
   ; of image. We calculate these in UV space. Note that these values are in the
   ; center of the pixel. We have to move them to the edge of the pixel below.
   ; The U direction corresponds to longitude, and the V direction to latitude.
   topv =   (uv[1,1] + uv[1,2]) * 0.5
   botv =   (uv[1,0] + uv[1,3]) * 0.5
   leftu =  (uv[0,0] + uv[0,1]) * 0.5
   rightu = (uv[0,2] + uv[0,3]) * 0.5

   ; Calculate the scales in UV directions. The variable s contains the image dimensions.
   xscale = (rightu-leftu) / (s[0]-1)
   yscale = (topv-botv) / (s[1]-1)

   ; Get the range of the UV box that envelopes the image.
   xrange = [uv[0,1] - (0.5 * xscale), uv[0,2] + (0.5 * xscale)]
   yrange = [uv[1,0] - (0.5 * yscale), uv[1,1] + (0.5 * yscale)]

   ; Decide on a position of the image in the window.
   pos = [0.05, 0.25, 0.95, 0.95]
   
   ; Update the map coordinate object.
   mapCoord -> SetProperty, XRANGE=xrange, YRANGE=yrange, POSITION=pos

   ; Display the image. The variable POS will change (probably) to keep the aspect.
   cgDisplay, 500, 500
   cgImage, scaledImage, POSITION=pos, /KEEP_ASPECT, /ERASE

   ; Set up a map coordinate space on top of the image in UV coordinates.
   mapCoord -> Draw
;
;   ; Add continent and state outlines, along with grid labels.
   cgMap_Continents, /HIRES, Color='medium gray', /USA, MAP=mapCoord
   cgMap_Grid, LATS=Indgen(8)*5+Round(MIN(latitude)), /LABEL, $
             LONS = Indgen(8)*10+Round(Min(longitude)), $
             COLOR='gray', MAP=mapCoord

   ; Add a colorbar. Non-linear scaling requires use of tick formatting function.
   cgColorbar, NColors=13, Bottom=1, Position=[pos[0], 0.1, pos[2], 0.15], $
      Divisions=14, Title='24 Hour Precipitation (mm)', AnnotateColor='black', $
      /Discrete, OOB_High=14, XTickFormat='PrecipMap_Annotate'
  
   ; Need output keywords?
   IF Arg_Present(data) THEN data = scaledImage
   IF Arg_Present(palette) THEN BEGIN
      TVLCT, r, g, b, /Get
      palette = [[r],[g],[b]]
   ENDIF
   IF Arg_Present(map) THEN map = mapCoord
        
END ;----------------------------------------------------------------------