PRO TiffOverlay

   filename = 'AF03sep15b.n16-VIg.tif'
   image = Read_Tiff(filename, GEOTIFF=geotag) ; GeoTIFF info in "geotag" structure.
   image = Reverse(image, 2) ; Reverse Y direction for display in IDL.

   ; Find the image dimensions. Will need later.
   s = Size(image, /Dimensions)

    ; Calculate corner points from GeoTIFF structure obtained from file.
   xscale  = geotag.ModelPixelScaleTag[0]
   yscale  = geotag.ModelPixelScaleTag[1]
   tp      = geotag.ModelTiePointTag[3:4]

   ; Origin in upper-left in TIFF file. I want to move to bottom left
   ; for my own sanity.
   xOrigin = tp[0]
   yOrigin = tp[1] - (yscale * s[1])
   xEnd = xOrigin + (xscale * s[0])
   yEnd = tp[1]
   
   ; Print the corner coordinates. They correspond exactly to what is in TIFF file
   ; according to listgeo from www.remotesensing.org.
   Print, 'Corner Coordinates:'
   Print, 'Upper Left:  (' + String(xOrigin, Format='(D0.3)') + $
          ', ' + String(yEnd, Format='(D0.3)') + ')'
   Print, 'Lower Left:  (' + String(xOrigin, Format='(D0.3)') + $
          ', ' + String(yOrigin, Format='(D0.3)') + ')'
   Print, 'Upper Right: (' + String(xEnd, Format='(D0.3)') + ', ' + $
           String(yEnd, Format='(D0.3)') + ')'
   Print, 'Lower Right: (' + String(xEnd, Format='(D0.3)') + ', ' + $
           String(yOrigin, Format='(D0.3)') + ')'
 
    ; Map the corners in lat/lon space.
    alberMap = MAP_PROJ_INIT('Albers Equal Area Conic', $
       DATUM=8, $ ; WGS84
       CENTER_LAT=geotag.PROJNATORIGINLATGEOKEY, $
       CENTER_LON=geotag.PROJNATORIGINLONGGEOKEY, $
       STANDARD_PAR1=geotag.PROJSTDPARALLEL1GEOKEY, $
       STANDARD_PAR2=geotag.PROJSTDPARALLEL2GEOKEY)
    lonlat = MAP_PROJ_INVERSE([xOrigin, xOrigin, xEnd, xEnd   ], $
                              [yOrigin, yEnd,    yEnd, yOrigin], $
                               MAP_STRUCTURE=alberMap)    
                                
    Print, ''
    Print, 'Image Corners in Latitude and Longitude'
    Print, 'Latitude:  ', Reform(lonlat[1,*])
    Print, 'Longitude: ', Reform(lonlat[0,*])
    Print, ''
    
     ; The limit is the min and max of these latitudes and longitudes.
     ; Alas, this doesn't give us the correct limit. We need an 8-element
     ; array to do so, rather than a 4-element array.
     minLat = Min(lonlat[1,*], MAX=maxLat)
     minLon = Min(lonlat[0,*], MAX=maxLon)
     limit = [minLat, minLon, maxLat, maxLon]

   ; Set up a map projection space, using the limit from above.
   ; The UV box that results will be slightly wrong.
   alberMap = MAP_PROJ_INIT('Albers Equal Area Conic', $     
      DATUM=8, $ ; WGS84
      CENTER_LAT=geotag.PROJNATORIGINLATGEOKEY, $
      CENTER_LON=geotag.PROJNATORIGINLONGGEOKEY, $
      STANDARD_PAR1=geotag.PROJSTDPARALLEL1GEOKEY, $
      STANDARD_PAR2=geotag.PROJSTDPARALLEL2GEOKEY, $
      LIMIT=limit)
      
   ; Erase the display and set up an image position.
   IF (!D.Flags AND 256) NE 0 THEN Erase, COLOR=FSC_COLOR('ivory')
   pos = [0,0,1,1.]

   ; Set up color table.
   LoadCT, 0, /Silent
   TVLCT, FSC_Color('ivory', /Triple), 0
   LoadCT, 33, NColors=253, Bottom=1, /Silent

   ; Scale the image.
   scaled = BytScl(image, Top=253)+1
   index = Where(scaled EQ 1)
   scaled[index] = 0B

   ; Display the image.
   Window, XSIZE=500, YSIZE=500, Title='Outline using UV_BOX with Incorrect Limits.', $
      XPOS=10, YPOS=10
   TVIMAGE, scaled, POSITION=pos, /NOINTERP, /KEEP_ASPECT
   
   ; Use the (slight wrong) UV_BOX of map structure to set up the data coordinate space.
   uv_box = alberMap.uv_box
   Print, 'MAP_PROJ_INIT UV_BOX:  ', uv_box
   Plot, uv_box[[0, 2]], uv_box[[1, 3]], Position=pos, $
          /Nodata, XStyle=5, YStyle=5, /NoErase
          
    ; Draw political boundaries and coasts.
    MAP_CONTINENTS,/HIRES, Color=FSC_Color('Dark Slate Blue'), Map_Structure=alberMap

    ; Add gridlines
    MAP_GRID,/LABEL, Color=FSC_Color('Slate Blue'), Map_Structure=alberMap
    
    ; For completeness, let's try to calculate the correct limit for the plot.
    lefty =  Abs(lonlat[1,0]-lonlat[1,1])/2.0 + Min([lonlat[1,0],lonlat[1,1]])
    leftx =  Abs(lonlat[0,0]-lonlat[0,1])/2.0 + Min([lonlat[0,0],lonlat[0,1]])
    topy =   Abs(lonlat[1,1]-lonlat[1,2])/2.0 + Min([lonlat[1,1],lonlat[1,2]])
    topx =   Abs(lonlat[0,1]-lonlat[0,2])/2.0 + Min([lonlat[0,1],lonlat[0,2]])
    righty = Abs(lonlat[1,2]-lonlat[1,3])/2.0 + Min([lonlat[1,2],lonlat[1,3]])
    rightx = Abs(lonlat[0,2]-lonlat[0,3])/2.0 + Min([lonlat[0,2],lonlat[0,3]])
    boty =   Abs(lonlat[1,3]-lonlat[1,0])/2.0 + Min([lonlat[1,3],lonlat[1,0]])
    botx =   Abs(lonlat[0,3]-lonlat[0,0])/2.0 + Min([lonlat[0,3],lonlat[0,0]])
    limit = [lefty, leftx, topy, topx, righty, rightx, boty, botx]

    alberMap = MAP_PROJ_INIT('Albers Equal Area Conic', $     
      DATUM=8, $ ; WGS84
      CENTER_LAT=geotag.PROJNATORIGINLATGEOKEY, $
      CENTER_LON=geotag.PROJNATORIGINLONGGEOKEY, $
      STANDARD_PAR1=geotag.PROJSTDPARALLEL1GEOKEY, $
      STANDARD_PAR2=geotag.PROJSTDPARALLEL2GEOKEY, $
      LIMIT=limit)
   
    ; Do do this correctly, we need to use the "uv_box" we have already calculated.
   Window, XSIZE=500, YSIZE=500, 4, Title='Attempt Using Correct 8-Element Limits', $
      XPOS=60, YPOS=60
   TVIMAGE, scaled, POSITION=pos, /NOINTERP, /KEEP_ASPECT

   ; We use the UV_BOX returned from MAP_PROJ_INIT with 8-element LIMIT.
   uv_box = alberMap.uv_box
   Print, '8-Element UV_BOX:      ', uv_box
   Plot, uv_box[[0, 2]], uv_box[[1, 3]], Position=pos, $
          /Nodata, XStyle=5, YStyle=5, /NoErase
          
    ; Draw political boundaries and coasts.
    MAP_CONTINENTS,/HIRES, Color=FSC_Color('Dark Slate Blue'), Map_Structure=alberMap

    ; Add gridlines
    MAP_GRID,/LABEL, Color=FSC_Color('Slate Blue'), Map_Structure=alberMap

    
    ; To get correct results, we need to use the "uv_box" we have already calculated.
   Window, XSIZE=500, YSIZE=500, 1, Title='Outline using Correct Calculated Limits.', $
      XPOS=120, YPOS=120
   TVIMAGE, scaled, POSITION=pos, /NOINTERP, /KEEP_ASPECT

   ; Use the calculated UV_BOX to set up the data coordinate space.
   uv_box = [xOrigin, yOrigin, xEnd, yEnd]
   Print, 'Calculated UV_BOX:     ', uv_box
   Plot, uv_box[[0, 2]], uv_box[[1, 3]], Position=pos, $
          /Nodata, XStyle=5, YStyle=5, /NoErase
          
    ; Draw political boundaries and coasts.
    MAP_CONTINENTS,/HIRES, Color=FSC_Color('Dark Slate Blue'), Map_Structure=alberMap

    ; Add gridlines
    MAP_GRID,/LABEL, Color=FSC_Color('Slate Blue'), Map_Structure=alberMap
    
    

   ; We can also get correct limits from MAP_PROJ_IMAGE.
   Window, XSIZE=500, YSIZE=500, 2, Title='Outline using UV_BOX of Map_Proj_Image', $
      XPOS=180, YPOS=180
   TVIMAGE, scaled, POSITION=pos, /NOINTERP, /KEEP_ASPECT

   ; Use the UV_BOX of Map_Proj_Image to set up the data coordinate space.
   warp = Map_Proj_Image(dist(10), $
      [xorigin, yorigin, xend, yend], $ ; image range in UV coords
      Image_Structure=alberMap, $
      Map_Structure=alberMap, $
      UVRANGE=uv_box)
   
   Print, 'MAP_PROJ_IMAGE UV_BOX: ', uv_box
   Plot, uv_box[[0, 2]], uv_box[[1, 3]], Position=pos, $
          /Nodata, XStyle=5, YStyle=5, /NoErase

    ; Put the world on top
    MAP_CONTINENTS,/HIRES, Color=FSC_Color('Dark Slate Blue'), Map_Structure=alberMap

    ; Add gridlines
    MAP_GRID,/LABEL, Color=FSC_Color('Slate Blue'), Map_Structure=alberMap
     
END