;+
; NAME:
;       READ_TOMS_AEROSOL
;
; PURPOSE:
;
;       This is a demonstration program, used in the TOMS Tutorial
;       (http://www.dfanning.com/graphics_tips/toms_tutorial.html),
;       to show how to write a program that looks the same on the
;       display and in a PostScript file. TOMS satellite data is
;       overlaid on a map projection.
;
; AUTHOR:
;
;       FANNING SOFTWARE CONSULTING
;       David Fanning, Ph.D.
;       1645 Sheely Drive
;       Fort Collins, CO 80526 USA
;       Phone: 970-221-0438
;       E-mail: davidf@dfanning.com
;       Coyote's Guide to IDL Programming: http://www.dfanning.com
;
; CATEGORY:
;
;       Graphics
;
; CALLING SEQUENCE:
;
;       Read_TOMS_Aerosol, filename
;
; ARGUMENTS:
;
;       filename:   The name of the TOMS data file. If not provided, the user
;                   will be prompted for the name of the file.
;
; KEYWORDS:
;
;       DATA:       Set this keyword to a named variable that on return will
;                   contain a 288 by 180 array of floating point values. Missing
;                   data is indicated by !Values.F_NAN.
;
;       HEADER:     Set this keyword to a named variable that on return will
;                   contain the first three lines of the TOMS data set.
;
;       LIMIT:      A four-element input vector that limits the extent of the map
;                   projection. Identical to the LIMIT keyword for MAP_SET.
;                   The form of the keyword is:
;
;                       LIMIT = [latmin, lonmin, latmax, lonmax]
;
;       PS:         Set this keyword to 1 if you wish to create PostScript output
;                   instead of sending graphics output to the display.
;
; COMMON BLOCKS:
;
;       None
;
; RESTRICTIONS:
;
;       TOMS data can be freely obtained from NASA at http://toms.gsfc.nasa.gov/.
;
;       Requires numerious files from the Coyote Library:
;
;           http://www.dfanning.com/programs/coyoteprograms.zip
;
; MODIFICATION HISTORY:
;
;       Written by David W. Fanning, 18 January 2006
;-
;
;###########################################################################
;
; LICENSE
;
; This software is OSI Certified Open Source Software.
; OSI Certified is a certification mark of the Open Source Initiative.
;
; Copyright © 2006 Fanning Software Consulting.
;
; This software is provided "as-is", without any express or
; implied warranty. In no event will the authors be held liable
; for any damages arising from the use of this software.
;
; Permission is granted to anyone to use this software for any
; purpose, including commercial applications, and to alter it and
; redistribute it freely, subject to the following restrictions:
;
; 1. The origin of this software must not be misrepresented; you must
;    not claim you wrote the original software. If you use this software
;    in a product, an acknowledgment in the product documentation
;    would be appreciated, but is not required.
;
; 2. Altered source versions must be plainly marked as such, and must
;    not be misrepresented as being the original software.
;
; 3. This notice may not be removed or altered from any source distribution.
;
; For more information on Open Source Software, visit the Open Source
; web site: http://www.opensource.org.
;
;###########################################################################
PRO Read_TOMS_Aerosol, filename, DATA=data, HEADER=header, LIMIT=limit, PS=ps

   ; Error handling
   Catch, theError
   IF theError NE 0 THEN BEGIN
      Catch, /Cancel
      void = Error_Message()
      IF N_Elements(lun) NE 0 THEN Free_Lun, lun ; Free file logical unit number.
      IF N_Elements(theState) NE 0 THEN Device, Decomposed=theState ; Preserve color mode.
      RETURN
   ENDIF

   ; Open the file for reading.
   IF N_Elements(filename) EQ 0 THEN BEGIN
      filename = Dialog_Pickfile(Filter='*.txt', Title='Select EP/TOMS file for READING...')
      IF filename EQ "" THEN RETURN
   ENDIF
   OpenR, lun, filename, /Get_Lun

   ; First three lines are header information.
   header = StrArr(3)
   READF, lun, header

   ; TOMS data is a floating point array of [288,180]. Data is on 1 degree
   ; grid in lat and 1.25 degree grid in lon. Data stored as 16-bit integers,
   ; and must be converted to floats by dividing by 10.0. NAN stored as 999.
   data = IntArr(288, 180)
   line = IntArr(288)
   FOR j=0,179 DO BEGIN
      READF, lun, line, Format='(1x,25I3)'
      data[0,j] = line
   ENDFOR
   Free_Lun, lun

   ; Find 999.
   indices = Where(data EQ 999, count)

   ; Convert data to floating point values.
   data = data / 10.0
   IF count GT 0 THEN data[indices] = !Values.F_NAN

   ; Create corresponding lat/lon vectors.
   lat = Scale_Vector(Findgen(180), -89.5, 89.5)
   lon = Scale_Vector(Findgen(288), -179.375, 179.375)

   ; If LIMITS are set, then find the proper subscripted values.
   IF N_Elements(limit) EQ 0 THEN BEGIN
      latsubs = [0,179]
      lonsubs = [0,287]
   ENDIF ELSE BEGIN
      IF N_Elements(limit) NE 4 THEN $
         Message, 'Usage: LIMIT=[latmin, lonmin, latmax, lonmax].'
      latsubs = Value_Locate(lat, [limit[0],limit[2]])
      lonsubs = Value_Locate(lon, [limit[1],limit[3]])
   ENDELSE

   ; Subscript the data
   lat_sub = lat[latsubs[0]:latsubs[1]]
   lon_sub = lon[lonsubs[0]:lonsubs[1]]
   data_sub = data[lonsubs[0]:lonsubs[1], latsubs[0]:latsubs[1]]

   ; Scale the data between 1 and 4.5 into 6 colors. Set the minimum value to 1.
   img = BytScl(data_sub, Min=1, Max=4.5, Top=6, /NAN) + 1B

   ; Locate the NANs in the subsetted data. These values will be set to background.
   indices = Where(Finite(data_sub) EQ 0, count)
   IF count GT 0 THEN img[indices] = 10

   ; Locate the values in the data that are less than 1. These values will be set to background.
   lessThanOne = Where(data_sub LT 1.0, ltcount)
   IF ltcount GT 0 THEN img[lessThanOne] = 10

   ; Device independent graphics is most easily accomplished with decomposed color.
   Device, Decomposed=0, Get_Decomposed=theState

   ; Load image display colors.
   LoadCT, 0, /Silent

   ; Aerosol colors, indices 1-7;
   TVLCT, FSC_Color(['LIGHT GRAY', 'DARK GRAY', 'YELLOW', 'GREEN YELLOW', $
                     'RED', 'MAROON', 'DARK RED'], /Triple), 1
   TVLCT, FSC_Color('BLACK', /Triple), 8 ; Drawing color.
   TVLCT, FSC_Color('SKY BLUE', /Triple), 250 ; Ocean color.
   TVLCT, FSC_Color('FOREST GREEN', /Triple), 251 ; Landmass color.
   TVLCT, FSC_Color('WHITE', /Triple), 252 ; Background color.

   ; The next part of the program takes place in the Z-graphics buffer.
   ; It is here we need to create a seamask and landmask image in order
   ; to create the proper background image for our graphics display.
   x_winsize = 750 ; This will turn int 7.5 inches in PostScript.
   y_winsize = 600 ; This will turn into 6.0 inches in PostScript.
   thisDevice = !D.Name
   Set_Plot, 'Z'
   Erase
   Device, Set_Resolution=[x_winsize, y_winsize], Z_Buffer=0

   ; Set up the map projection center coordinates.
   lat_center = (Max(lat_sub) - Min(lat_sub)) / 2 + Min(lat_sub)
   lon_center = (Max(lon_sub) - Min(lon_sub)) / 2 + Min(lon_sub)

   ; Set up the map projection coordinate space, and draw the continents in a filled color.
   Map_Set, /Cylindrical, lat_center, lon_center, Position=[0.1, 0.25, 0.9, 0.9], $
      Limit=[Min(lat_sub), Min(lon_sub), Max(lat_sub), Max(lon_sub)], /NoBorder, /Isotropic
   Map_Continents, /Fill, Color=251

   ; Warp the image.
   warp = Map_Image(img, xx, yy, xs, ys, LatMin=Min(lat_sub), LonMin=Min(lon_sub), $
      LatMax=Max(lat_sub), LonMax=Max(lon_sub), Compress=1, Missing=252)

   ; Determine the size of the warped image.
   dims = Size(warp, /Dimensions)

   ; Grab this portion of the window to obtain sea and land mask images.
   snap = TVRD(xx, yy, dims[0], dims[1])
   seamask = snap EQ 0
   landmask = snap NE 0
   seacolor = 250
   seaimage = seamask * (BytArr(dims[0],dims[1]) + seacolor)
   landcolor = 251
   landimage = landmask * (BytArr(dims[0],dims[1]) + landcolor)

   ; Create the background image.
   backgroundImage = seaimage > landimage

   ; Back to the current device.
   Set_Plot, thisDevice

   ; Prepare display image. Set all pixels in the warped image with value 0
   ; are set to the background image color.
   indices = Where(warp EQ 10, count)
   IF count GT 0 THEN warp[indices] = backgroundImage[indices]


   ; Set up the appropriate window size for PostScript or display.
   IF Keyword_Set(PS) THEN BEGIN
      ps_xsize = x_winsize/100
      ps_ysize = y_winsize/100
      keywords = PSConfig(Cancel=cancelled, XSize=ps_xsize, YSize=ps_ysize, $
         Inches=1, XOffset=(8.5 - ps_xsize)/2.0, YOffset=(11.0 - ps_ysize)/2.0, $
         Filename='aerosol.ps')

      ; If the user cancelled, then forget about it.
      IF cancelled THEN BEGIN
         Device, Decomposed=theState
         RETURN
      ENDIF
      thisDevice = !D.Name
      Set_Plot, 'PS', /Copy
      Device, _Extra=keywords

      ; Have to perform the image warp again in order to get the proper size and start information.
      void = Map_Image(img, xx, yy, xs, ys, LatMin=Min(lat_sub), LonMin=Min(lon_sub), $
         LatMax=Max(lat_sub), LonMax=Max(lon_sub))

   ENDIF ELSE BEGIN

      Window, /Free, Title='TOMS Aerosol Data Plot', XSize=x_winsize, YSize=y_winsize
      Erase, Color=252

   ENDELSE

   ; Display the image.
   TV, warp, xx, yy, XSize=xs, YSize=ys

   ; Add continents, countries, and grids.
   Map_Continents, Color=8
   Map_Continents, Color=8, Countries=1
   Map_Grid, Color=8

   ; Add a colorbar.
   xstart = Float(xx)/!D.X_Size
   xend = Float(xs)/!D.X_Size + xstart
   yend = Float(yy)/!D.Y_Size -  0.10
   ystart = yend - 0.075
   tickvalues = ['1.0', '1.5', '2.0', '2.5', '3.0', '3.5', '4.0', '4.5+']
   Colorbar, NColors=7, Bottom=1, Range=[1.0, 4.5], Divisions=7, Format='(F3.1)', $
      Position=[xstart, ystart, xend, yend], XMinor=1, XTicklen=1, $
      Title='Aerosol Index', Font=0, Color=8, Ticknames=tickvalues

   ; Finish.
   IF Keyword_Set(PS) THEN BEGIN
      Device, /Close_File
      Set_Plot, thisDevice
   ENDIF

   ; Restore system state.
   Device, Decomposed=theState

END