PRO GSHHS_PLOT, FILE, LEVEL=LEVEL, AREA=AREA, COLOR=COLOR, $
  BACKGROUND=BACKGROUND, FILL=FILL, MAP_STRUCTURE=MAP_STRUCTURE, _EXTRA=_EXTRA

;+
; DESCRIPTION:
;    Plot coastlines from a GSHHS coastline file.
;
;    GSHHS stands for the Global Self-consistent, Hierarchical,
;    High-resolution Shoreline Database assembled by Paul Wessel
;    and Walter Smith. For more information:
;    http://www.ngdc.noaa.gov/mgg/shorelines/gshhs.html
;
;    The binary data files may be downloaded from the site above:
;    gshhs_f.b	Full resolution data
;    gshhs_h.b	High resolution data
;    gshhs_i.b	Intermediate resolution data
;    gshhs_l.b	Low resolution data
;    gshhs_c.b	Crude resolution data
;
;    NOTE: These are binary files in big-endian byte order.
;    Bytes will be flipped automatically if this program is run on a little-endian platform.
;
; USAGE:
;    GSHHS_PLOT, FILE
;
; INPUT PARAMETERS:
;    FILE          Full path and name of the GSHHS coastline file.
;
; OUTPUT PARAMETERS:
;    None.
;
; OPTIONAL KEYWORDS:
;    LEVEL         The level of geographical features plotted (default = 1).
;                  1 = land
;                  2 = lake
;                  3 = island in lake
;                  4 = pond in island in lake
;                  The selected level and all lower levels are plotted.
;    AREA          Minimum area of features to be plotted (sq. km; default = 500).
;    COLOR         Color table index for plotting (default = !D.TABLE_SIZE - 1).
;    BACKGROUND    Color table index for lakes or ponds when FILL is active (default = 0).
;    FILL          If set, land masses are filled with COLOR, and inland water
;                  bodies are filled with BACKGROUND.
;    MAP_STRUCTURE Set this keyword to a !MAP structure, as returned from MAP_PROJ_INIT. 
;
;    Also accepts any keywords accepted by PLOTS or POLYFILL.
;
; MODIFICATION HISTORY:
; Liam.Gumley@ssec.wisc.edu
; Copyright (C) 2006 Liam E. Gumley
;
; This program is free software; you can redistribute it and/or
; modify it under the terms of the GNU General Public License
; as published by the Free Software Foundation; either version 2
; of the License, or (at your option) any later version.
;
; This program is distributed in the hope that it will be useful,
; but WITHOUT ANY WARRANTY; without even the implied warranty of
; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
; GNU General Public License for more details.
;
; You should have received a copy of the GNU General Public License
; along with this program; if not, write to the Free Software
; Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
;-

;- Check arguments
if (n_elements(file) eq 0) then message, 'Argument FILE is undefined'
if (n_elements(level) eq 0) then level = 1
if (n_elements(area) eq 0) then area = 500.0
if (n_elements(color) eq 0) then color = !d.table_size - 1
if (n_elements(background) eq 0) then background = 0

;- Open the file (note: GSHSS binary data files are big-endian)
openr, lun, file, /get_lun, /swap_if_little_endian

;- Define GSHHS polygon header (v1.3 format)
header = {$
  id:0L, $
  npoints:0L, $
  level:0L, $
  west:0L, $
  east:0L, $
  south:0L, $
  north:0L, $
  area:0L, $
  version:0L, $
  greenwich:0S, $
  source:0S}
          
;- Read each polygon
while (eof(lun) ne 1) do begin

  ;- Read this polygon header
  readu, lun, header
  
  ;- Read this polygon record
  data = lonarr(2, header.npoints, /nozero)
  readu, lun, data

  ;- Check level and area of this polygon
  this_level = header.level
  this_area = header.area * 0.1
  if (this_level le level) and (this_area ge area) then begin

    ;- Extract longitude and latitude data
    lon = data[0, *] * 1.0e-6
    lat = data[1, *] * 1.0e-6

    ;- If MAP_PROJ structure is defined, convert to data coordinates
    if (n_elements(map_structure) ne 0) then begin
      xy = map_proj_forward(lon, lat, map_structure=map_structure)
      x = xy[0, *]
      y = xy[1, *]
    endif else begin
      x = lon
      y = lat
    endelse

    ;- Check if this polygon should be displayed
    result = convert_coord(x, y, /data, /to_normal)
    xnorm = result[0, *]
    ynorm = result[1, *]
    loc = (xnorm ge 0.0 and xnorm le 1.0) and (ynorm ge 0.0 and ynorm le 1.0)
    flag = max(loc)
    
    ;- Plot the polygon
    if (flag) then begin
      if (keyword_set(fill) eq 0) then begin
        plots, x, y, color=color, _extra=_extra
      endif else begin
        if (header.level eq 1) or (header.level eq 3) then begin
          polyfill, x, y, color=color, _extra=_extra
        endif else begin
          polyfill, x, y, color=background, _extra=_extra
        endelse    
      endelse
    endif
    
  endif
      
endwhile

;- Close the file
free_lun, lun

END