PRO GOES_EXAMPLE, filename

    ; Restore GOES example data. Original image and navigation data obtained 
    ; from http://goes.gsfc.nasa.gov/goeseast-lzw/peru/. Navigation arrays
    ; saved a long integers. Divide by 10000 to get latitude and longitude
    ; floats. Variables: peruimage, peru_lat, peru_lon.
    
    ; Need a save file?
    IF N_Elements(filename) EQ 0 THEN BEGIN
        filename = Dialog_Pickfile(FILTER='*.sav', TITLE='Open GOES Example Data...')
        IF filename EQ "" THEN RETURN
    ENDIF
    Restore, filename
          
    ; Restore the data to floats. 
    peru_lat = Temporary(peru_lat) / 10000.
    peru_lon = Temporary(peru_lon) / 10000.
    Help, peruimage, peru_lat, peru_lon
    
    ; Get the image dimensions for subsequent processing.
    s = Size(peruimage, /DIMENSIONS)
    Print, 'Dimensions of Original Image: ', s
    
    ; Locate the longitude and latitude in the center of the image.
    ; We will center the map projection here.
    center_lat = peru_lat[s[0]/2, s[1]/2]
    center_lon = peru_lon[s[0]/2, s[1]/2]
    Print, 'Center Latitude Original Image: ', center_lat
    Print, 'Center Longitude Original Image:', center_lon
      
    ; Load colors.
    LoadCT, 0, /Silent
        
    ; Open a display window.
    Window, XSIZE=800, YSIZE=575, TITLE='GOES Image in Albers Projection', $
        XPOS=20, YPOS=20
    XYOUTS, 0.5, 0.5, /NORMAL, FONT=0, ALIGN=0.5, $
        'The image is being warped. Please do not close this window.'
        
    ; Find the lat/lon values at the sides of the image, so you
    ; can specify an 8-element LIMIT in MAP_SET.
    lon_cindex = s[0]/2
    lat_cindex = s[1]/2
    
    ; Side of original image.
    left_lon = peru_lon[0, lat_cindex]
    top_lon = peru_lon[lon_cindex, s[1]-1]
    right_lon = peru_lon[s[0]-1, lat_cindex]
    bottom_lon = peru_lon[lon_cindex, 0]
    
    left_lat = peru_lat[0, lat_cindex]
    top_lat = peru_lat[lon_cindex, s[1]-1]
    right_lat = peru_lat[s[0]-1, lat_cindex]
    bottom_lat = peru_lat[lon_cindex, 0]
    
    Print, 'Sides of original image in lat/lon space.'
    Print, Reform([[left_lon, top_lon, right_lon, bottom_lon], $
                   [left_lat, top_lat, right_lat, bottom_lat]], 2, 4)
    
    ; Corners of the original image.
    ulc_lon = peru_lon[0,s[1]-1]
    urc_lon = peru_lon[s[0]-1, s[1]-1]
    lrc_lon = peru_lon[s[0]-1, 0]
    llc_lon = peru_lon[0,0]
    
    ulc_lat = peru_lat[0,s[1]-1]
    urc_lat = peru_lat[s[0]-1, s[1]-1]
    lrc_lat = peru_lat[s[0]-1, 0]
    llc_lat = peru_lat[0,0]


    Print, 'Corners of original image in lat/lon space.'
    corners = [[ulc_lon, urc_lon, lrc_lon, llc_lon], $
               [ulc_lat, urc_lat, lrc_lat, llc_lat]]
    Print, Transpose(corners)
                   
    ; Create an Albers Equal Area Conic map projection space.
    Map_Set, center_lat, center_lon, /Albers, /NOBORDER, $
        XMARGIN=0, YMARGIN=0, STANDARD_PARALLEL=[-19, 21], $
        LIMIT=[left_lat, left_lon, top_lat, top_lon, $
               right_lat, right_lon, bottom_lat, bottom_lon], /NOERASE
               
    ; Warp the image and the navigation arrays.
    warp = Map_Patch(peruimage, peru_lon, peru_lat)
    warp_lat = Map_Patch(peru_lat, peru_lon, peru_lat)
    warp_lon = Map_Patch(peru_lon, peru_lon, peru_lat)
               
    ; The warping process can add strangeness to the pixels on
    ; the outside edge of the arrays. Lots of zeros, etc. To avoid
    ; problems, we just throw away the outer pixels.
    s = Size(warp, /DIMENSIONS)
    warp =         warp[1:s[0]-2, 1:s[1]-2]
    warp_lat = warp_lat[1:s[0]-2, 1:s[1]-2] 
    warp_lon = warp_lon[1:s[0]-2, 1:s[1]-2]  
    center_lat = warp_lat[s[0]/2, s[1]/2]
    center_lon = warp_lon[s[0]/2, s[1]/2]
    
    Print, 'Center of warped image: ', center_lon, center_lat
    
    
    s = Size(warp, /DIMENSIONS)
    Print, 'Dimensions of warped image: ', s
    
    ; Side of warped image.
    left_lon = warp_lon[0, lat_cindex]
    top_lon = warp_lon[lon_cindex, s[1]-1]
    right_lon = warp_lon[s[0]-1, lat_cindex]
    bottom_lon = warp_lon[lon_cindex, 0]
    
    left_lat = warp_lat[0, lat_cindex]
    top_lat = warp_lat[lon_cindex, s[1]-1]
    right_lat = warp_lat[s[0]-1, lat_cindex]
    bottom_lat = warp_lat[lon_cindex, 0]

    Print, 'Sides of warped image in lat/lon space.'
    Print, Reform([[left_lon, top_lon, right_lon, bottom_lon], $
                   [left_lat, top_lat, right_lat, bottom_lat]], 2, 4)

        ; Corners of the warped image.
    ulc_lon = warp_lon[0,s[1]-1]
    urc_lon = warp_lon[s[0]-1, s[1]-1]
    lrc_lon = warp_lon[s[0]-1, 0]
    llc_lon = warp_lon[0,0]
    
    ulc_lat = warp_lat[0,s[1]-1]
    urc_lat = warp_lat[s[0]-1, s[1]-1]
    lrc_lat = warp_lat[s[0]-1, 0]
    llc_lat = warp_lat[0,0]
 
    Print, 'Corners of warped image in lat/lon space.'
    corners = [[ulc_lon, urc_lon, lrc_lon, llc_lon], $
               [ulc_lat, urc_lat, lrc_lat, llc_lat]]
    Print, Transpose(corners)
    
    ; Now create the map projection space for the warped image, using the
    ; new parameters from above.
    Map_Set, center_lat, center_lon, /Albers, /NOBORDER, $
        XMARGIN=0, YMARGIN=0, STANDARD_PARALLEL=[-19, 21], $
        LIMIT=[left_lat, left_lon, top_lat, top_lon, $
               right_lat, right_lon, bottom_lat, bottom_lon], /NOERASE
    
    ; Display the warped image in the window.
    TVIMAGE, warp, /NOINTERP
    
    ; Add annotation to the image. Be sure to use high resolution maps.
    Map_Continents, COLOR=FSC_COLOR('Pale Goldenrod'), /HIRES
    Map_Continents, /COUNTRIES, COLOR=FSC_COLOR('Pale Goldenrod'), /HIRES
    Map_Grid, /LABEL, COLOR=FSC_COLOR('Light Gray')
 
    ; Now we wish to save the warped and map projected image as a GEOTIFF file.
    ; We have to specify scales and tie points in the UV coordinate space. So
    ; we have to project information from lat/lon space to UV space (a forward
    ; transformation. Use a GCTP map projection to do so.
     alberMap = MAP_PROJ_INIT('Albers Equal Area Conic', $
       DATUM=8, $ ; WGS84
       CENTER_LAT=center_lat, $
       CENTER_LON=center_lon, $
       STANDARD_PAR1=-19, $
       STANDARD_PAR2=20)
 
 
    ; We want to transform the sides of the image from lat/lon to UV space.
    uv = MAP_PROJ_FORWARD([left_lon, top_lon, right_lon, bottom_lon], $
                          [left_lat, top_lat, right_lat, bottom_lat], $
                          MAP_STRUCTURE=alberMap)
                          
    Print, 'Sides of warped image in UV coordinates'
    Print, uv
    
  
    ; Calculate the tie point (upper left corner) and scale
    ; factors in UV space.
    s = Size(warp, /DIMENSIONS)
    xscale = Abs(uv[0,0] - uv[0,2]) / (s[0] )
    yscale = Abs(uv[1,1] - uv[1,3]) / (s[1])
    tp = [uv[0,0], uv[1,1]]
;    xscale = Abs(uvc[0,0] - uvc[0,1]) / (s[0])
;    yscale = Abs(uvc[1,1] - uvc[1,2]) / (s[1])
;    tp = [uvc[0,0], uvc[1,0]]
    Print, 'Scales For GeoTiff File:    ', xscale, yscale
    Print, 'Tie Point for GeoTiff File: ', tp
  
    ; Fill out the geotag structure with the proper information.
    ; See the GEOTIFF web page for definitions: 
    ; http://www.remotesensing.org/geotiff/geotiff.html
    print, 'Centers:', center_lon, center_lat
    geotag = { MODELPIXELSCALETAG: [xscale, yscale, 0], $
               MODELTIEPOINTTAG: [ 0, 0, 0, tp, 0], $
               GTMODELTYPEGEOKEY: 1, $
               GTRASTERTYPEGEOKEY: 1, $
               GEOGRAPHICTYPEGEOKEY: 4326, $
               GEOGLINEARUNITSGEOKEY: 9001, $
               GEOGANGULARUNITSGEOKEY: 9102, $
               PROJECTEDCSTYPEGEOKEY: 32767, $
               PROJECTIONGEOKEY: 32767, $
               PROJCOORDTRANSGEOKEY: 11, $
               PROJLINEARUNITSGEOKEY: 9001, $
               PROJSTDPARALLEL1GEOKEY: -19.00000, $
               PROJSTDPARALLEL2GEOKEY: 21.000000, $
               PROJNATORIGINLONGGEOKEY: center_lon, $
               PROJNATORIGINLATGEOKEY: center_lat, $
               PROJFALSEEASTINGGEOKEY: 0, $
               PROJFALSENORTHINGGEOKEY: 0 }
               
    ; Write the GeoTiff file.
    Write_TIFF, 'goes_example.tif', Byte(Reverse(warp, 2)), GEOTIFF=geotag
       
;###################################################################################


    ; Open and read the GeoTiff image. Be sure to get the geotag structure.
    tiffimage = 'goes_example.tif'
    image = Read_Tiff(tiffimage, GEOTIFF=geotag)
    image = Reverse(image, 2)
    
   ; Find the image dimensions.
   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]
   Print, 'Scales in Display Program', xscale, yscale
   Print, 'Tie point in Display Program', tp

   ; Origin in upper-left in TIFF file. 
   xOrigin = tp[0]
   yOrigin = tp[1]
   xEnd = xOrigin + (xscale * s[0])
   yEnd = yOrigin - (yscale * s[1])
   
   Print, 'Image Corners in UV Space'
   corners = Fltarr(2,4)
   corners[0,*] = [xorigin, xend, xend, xorigin]
   corners[1,*] = [yorigin, yorigin, yend, yend] 
   print, corners, format='(2f14.1)'
   
   ; Map the corners in lat/lon space by inverse transformation.
    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, xEnd, xEnd, xOrigin], $
                              [yEnd, yEnd, yOrigin, yOrigin], $
                               MAP_STRUCTURE=alberMap) 
                               
    Print, 'Image Corners in Lat/Lon Space'
    Print, lonlat

    ; Create a window and display the image.
    s = Size(image, /Dimensions)
    Window, 1, XSIZE=s[0], YSIZE=s[1], TITLE='GeoTiff Output', XPOS=70, YPOS=70
    pos = [0, 0, 1, 1]
    TVIMAGE, image, POSITION=pos, /NOINTERP

    ; Use the image corners to set up the data coordinate space.
    Plot, [xOrigin, xEnd], [yEnd, yOrigin], Position=pos, $
          /Nodata, XStyle=5, YStyle=5, /NoErase
          
    ; Add map annotation.
    Map_Continents, COLOR=FSC_COLOR('Pale Goldenrod'), /HIRES, MAP_STRUCTURE=alberMap
    Map_Continents, /COUNTRIES, COLOR=FSC_COLOR('Pale Goldenrod'), /HIRES, MAP_STRUCTURE=alberMap
    Map_Grid, /LABEL, COLOR=FSC_COLOR('Light Gray'), MAP_STRUCTURE=alberMap, $
        LATS=Indgen(10) * 4 + (-20), LONS=Indgen(12) * 5 + (-105)

END