;---------------------------------------------------------------------------
;+
; NAME:
;       PIMAGE
; PURPOSE:
;       Display or plot a color image with axis and a color bar, using
;       one of several projections (default is cylindrical equidistant).
; LAST UPDATE: 11 July 2000 
; CATEGORY:
; CALLING SEQUENCE:
;       pimage      -> Displays help information for pimage.
;       pimage, a
; MODIFICATION HISTORY:
;       W. Berg - April 3, 1995, Original ce_image.pro code  
;       D. Anderson, Feb 1999, Added _extra structure to pass add'l args 
;       D. Anderson, Mar 1999, Added projection options 
;       D. Anderson, June 1999, Added GIF options 
;       D. Anderson, July 1999, Added region and (modified) smooth options 
;       D. Anderson, Aug 1999, Device available colors check 
;       D. Anderson, Aug 1999, Change label defaults and imgsize handling 
;       G. Wick, Apr 2000, Added /islands flag for hi-resolution map and coasts
;       D. Anderson, July 2000, Modified GIF handling, including !GFILE system var 
;       G. Dargaud, Jan 2001, Added PNG options
;
; VARIABLE DESCRIPTIONS:
;   dcnt   Count of valid data values in passed array 
;   dndex  Array index to valid data values in passed array 
;   px     Lower-left x coordinate for plotted image, in full page coordiates
;   py     Lower-left y coordinate for plotted image, in full page coordiates
;   pxc    Lower-left x coordinate for plotted colorbar, in full page coordiates
;   pyc    Lower-left y coordinate for plotted colorbar, in full page coordiates
;   xs     Width of plotted image, in full page coordiates
;   ys     Height of plotted image, in full page coordiates
;   xsc    Width of plotted colorbar, in full page coordiates
;   ysc    Height of plotted colorbar, in full page coordiates
;---------------------------------------------------------------------------
pro pimage, a, title=title, subtitle=sttl, units=units, project=prj, $ 
            xlabel=xlab, ylabel=ylab, vmin=vmn, vmax=vmx, region=region, $
            minlat=minlat, maxlat=maxlat, minlon=minlon, maxlon=maxlon, $
            nxminor=mxt, nyminor=myt, mis_val=misval, msk_val=mskval, $
            flg_val=flgval, mis_col=miscol, msk_col=mskcol, flg_col=flgcol, $
            multi=multi, $
            gifout=gifout, giffile=gfile, pngout=pngout, pngfile=pfile, $
            ps_port=psp, ps_land=psl, eps_port=epsp, eps_land=epsl, $ 
            coast=coast, usa=usa, ccolor=ccol, gcolor=gcol, cthick=cth, $ 
            gthick=gth, ctype=ctyp, gtype=gtyp, exp=exp, left_off=left_off, $
            right=right, top=top, bottom_off=bottom_off, cbar_off=cbar_off, $
            label=label, order=order, back_col=bcol, axis_col=acol, lsize=ls, $
            tsize=ts, position=pos, col_file=cfile, ps_file=psfile, $
            tick_off=tick_off, smooth=smooth, aspect=aspect, xtickname=xtn, $
            ytickname=ytn, ctickname=ctn, xtickv=xtv, ytickv=ytv, $
            ncolors=nc, vertical=vbar, horizontal=hbar, $
            imgpos=imgpos, islands=islands, _extra=xtra
 
if (n_params(0) EQ 0) then begin
  print,' PIMAGE   -> Displays help information for image.'
  print,' PIMAGE,a -> Displays a color image with axis and a color bar.'
  print,' Keywords:'
  print,'   ASPECT=aspect  Aspect ratio for the plot (def=(maxlon-minlon)/(maxlat-minlat).'
  print,'   SCALE=scale    Scale factor for the data.'
  print,'   OFFSET=offset  Offset factor for the data.'
  print,'   TITLE=title    Plot title'
  print,'   SUBTITLE=sttl  Secondary plot title'
  print,'   UNITS=units    Units label for the color bar'
  print,'   PROJECT=prj    Map projection (default: cylindrical equidistant)' 
  print,'                  cyl  Cylindrical Equidistant' 
  print,'                  mer  Mercator' 
  print,'                  mol  Mollweide' 
  print,'                  sat  Satellite'
  print,'                  polN Polar stereographic, north pole' 
  print,'                  polS Polar stereographic, south pole'  
  print,'                  xy   XY plot (no map)'
  print,'   SUBLAT=sublat  Subsatellite latitude for satellite projection' 
  print,'   SUBLNG=sublng  Subsatellite longitude for satellite projection' 
  print,'   SATPOS=satpos  Satellite p, omega, gamma' 
  print,'   XLABEL=xlab    X axis label.'
  print,'   YLABEL=ylab    Y axis label.'
  print,'   VMIN=vmn       Minimum data value.'
  print,'   VMAX=vmx       Maximum data value.'
  print,'   CMIN=cmn       Color value corresponding to the minimum data value.'
  print,'   CMAX=cmx       Color value corresponding to the maximum data value.'
  print,'   MINLAT=minlat  Minimum latitude value.'
  print,'   MAXLAT=maxlat  Maximum latitude value.'
  print,'   MINLON=minlon  Minimum longitude value.'
  print,'   MAXLON=maxlon  Maximum longitude value.'
  print,'   REGION=region  [minlat,maxlat,minlon,maxlon] array.' 
  print,'   LEVEL=nl       Number of the level to display (1..n).'
  print,'   TIME=nt        Number of the time band to display (1..n).'
  print,'   VAR=nv         Number of the variable to display (1..n).'
  print,'   NXTICK=nxt     Number of major tick mark intervals along x'
  print,'   NYTICK=nyt     Number of major tick mark intervals along y.'
  print,'   NCTICK=nct     Number of tick mark intervals for the color bar.'
  print,'   NXMINOR=mxt    Minor tick intervals per major tick along x.'
  print,'   NYMINOR=myt    Minor tick intervals per major tick along y.'
  print,'   XTICKNAME=xtn  Array of tick names for the x axis.'
  print,'   YTICKNAME=ytn  Array of tick names for the y axis.'
  print,'   CTICKNAME=ctn  Array of tick names for the colorbar axis.'
  print,'   XTICKV=xtv     Array of tick positions for the x axis.'
  print,'   YTICKV=ytv     Array of tick positions for the y axis.'
  print,'   CTICKLEN=ctl   Multiplier to adjust the color bar tick length.'
  print,'   XTICKLEN=xtl   Multiplier to adjust the xaxis bar tick length.'
  print,'   YTICKLEN=ytl   Multiplier to adjust the yaxis tick length.'
  print,'   TICKLEN=ctl    Multiplier to adjust the tick mark lengths.'
  print,'   TSIZE=ts       Character size for title.'
  print,'   LSIZE=ls       Character size for labels.'
  print,'   BACK_COL=bcol  Color for the background.'
  print,'   AXIS_COL=acol  Color for axis lines and labels.'
  print,'   MIS_VAL=misval The missing data value. (def=-999)'
  print,'   MSK_VAL=mskval The data value to mask out. (def=-998)'
  print,'   FLG_VAL=flgval The data value to flag. (def=-997)'
  print,'   MIS_COL=miscol The color value for missing data. (def=0)'
  print,'   MSK_COL=mskcol The color value for mask data. (def=2)'
  print,'   FLG_COL=flgcol The color value for flagged data. (def=4)'
  print,'   CCOLOR=ccol    Color for coastal outlines.'
  print,'   GCOLOR=gcol    Color for grid lines.' 
  print,'   CTHICK=cth     Thickness of coastal outlines.'
  print,'   GTHICK=gth     Thickness of grid lines.' 
  print,'   ATHICK=ath     Thickness of axis lines.' 
  print,'   CTYPE=ctyp     Style of coastal outlines.' 
  print,'   GTYPE=gtyp     Style of grid lines.' 
  print,'   POSITION=pos   Position of plot within window (x1,y1,x2,y2)'
  print,'                  (def=[0,0,1,1]->whole window).'
  print,'   IMGPOS=imgpos  Position of the image within the plot window.'
  print,'                  (image x1,y1,wid,ht & colorbar x1,y1,wid,ht)' 
  print,'   COL_FILE=cfile Name of the ascii color table file.'
  print,'   PS_FILE=psfile Name of the output PostScript file.'
  print,'   GIFFILE=gfile  Name of the output GIF file.'
  print,'   PNGFILE=pfile  Name of the output PNG file.'
  print,'   NPREC=nprec    Number of decimal points precision for lat/lon labels.'
  print,'   XWS=xws        Horizontal size of the output window in pixels.'
  print,'   YWS=yws        Vertical size of the output window in pixels.'
  print,'   SMOOTH=bcar    smooth image using boxcar average (1,3,5, etc.)'
  print,'   /MULTI         Set flag for doing multiple plots (do not set for last plot).'
  print,'   /PS_PORT       Write output to PostScript file in portrait mode.'
  print,'   /EPS_PORT      Write output to Encapsulated PS file, portrait mode.'
  print,'   /PS_LAND       Write output to PostScript file in landscape mode.'
  print,'   /EPS_LAND      Write output to Encapsulated PS file, landscape mode.'
  print,'   /GIFOUT        Write output as a GIF image.' 
  print,'   /PNGOUT        Write output as a PNG image.' 
  print,'   /COAST         Draw coastlines over the image.'
  print,'   /EXP           Exponential color scale (2^ncolors)'
  print,'   /BOTTOM_OFF    Turn off axis on the bottom of plot.'
  print,'   /LEFT_OFF      Turn off axis on the left side of plot.'
  print,'   /TOP           Add and label axis on top of plot.'
  print,'   /RIGHT         Add and label axis on right side of plot.'
  print,'   /TICK_OFF      Turn off the tick marks on the axis and color bar.'
  print,'   /CBAR_OFF      Turn off the color bar.'
  print,'   /VERTICAL      Force the color bar to be vertically oriented.'
  print,'   /HORIZONTAL    Force the color bar to be horizontally oriented.'
  print,'   /LABEL         Label axis E/W and N/S (default is  +/-)'
  print,'   /ORDER         Reverse order of image (top to bottom)'
  print,'   /USA           Draw U.S. state outlines.'
  print,'   /ISLANDS       Use high-resolution map and draw islands.'
  print,'   /VERBOSE       Print status information.'

  return
endif

; Default variable values

 cmn=5
 cmx=255
 nl=1 
 nt=1
 nv=1
 nct=5
 nxt=6
 nyt=6
 ctl=1.0
 xtl=1.0
 ytl=1.0
 offset=0.0
 sublat=0.0
 sublng=-135.0
 scale=1.0
 screenflag=0
 xws=800
 yws=600

; Check for above variables passed as "extra" (non-explicit) arguments  

 varcnt = n_tags(xtra)
 if varcnt gt 0 then begin 
   varstr=tag_names(xtra)
   for i = 1, varcnt do begin 
     if strpos(varstr(i-1), 'LEVEL') gt -1 then nl=xtra.level 
     if strpos(varstr(i-1), 'TIME') gt -1 then nt=xtra.time
     if strpos(varstr(i-1), 'VAR') gt -1 then nv=xtra.var 
     if strpos(varstr(i-1), 'NCTICK') gt -1 then nct=xtra.nctick 
     if strpos(varstr(i-1), 'NXTICK') gt -1 then nxt=xtra.nxtick 
     if strpos(varstr(i-1), 'NYTICK') gt -1 then nyt=xtra.nytick 
     if strpos(varstr(i-1), 'SUBLAT') gt -1 then sublat=xtra.sublat 
     if strpos(varstr(i-1), 'SUBLNG') gt -1 then sublng=xtra.sublng 
     if strpos(varstr(i-1), 'SCALE') gt -1 then scale=xtra.scale 
     if strpos(varstr(i-1), 'OFFSET') gt -1 then offset=xtra.offset 
     if strpos(varstr(i-1), 'CMIN') gt -1 then cmn=xtra.cmin 
     if strpos(varstr(i-1), 'CMAX') gt -1 then cmx=xtra.cmax 
     if strpos(varstr(i-1), 'XWS') gt -1 then xws=xtra.xws 
     if strpos(varstr(i-1), 'YWS') gt -1 then yws=xtra.yws 
     if strpos(varstr(i-1), 'CTICKLEN') gt -1 then ctl=xtra.cticklen
     if strpos(varstr(i-1), 'XTICKLEN') gt -1 then xtl=xtra.xticklen
     if strpos(varstr(i-1), 'YTICKLEN') gt -1 then ytl=xtra.yticklen
     if strpos(varstr(i-1), 'NPREC') gt -1 then nprec=xtra.nprec
     if strpos(varstr(i-1), 'XSIZE') gt -1 then xsz=xtra.xsize
     if strpos(varstr(i-1), 'YSIZE') gt -1 then ysz=xtra.ysize
     if strpos(varstr(i-1), 'ATHICK') gt -1 then ath=xtra.athick
     if strpos(varstr(i-1), 'VERBOSE') gt -1 then verbose=xtra.verbose
     if strpos(varstr(i-1), 'TLEFT') gt -1 then tleft=xtra.tleft
     if strpos(varstr(i-1), 'TAPER') gt -1 then taper=xtra.taper
     if strpos(varstr(i-1), 'TCOL') gt -1 then tcol=xtra.tcol
   end
 endif

; Read in default header values from common block

a2 = hdr_intf(a, title, sttl, units, minlat, maxlat, minlon, maxlon, misval, $
              mskval, order, nl, nt, nv, verbose=verbose)

;-------  Set remaining defaults  ---------------

 if n_elements(title) eq 0 then title=''
 if n_elements(sttl) eq 0 then sttl=''
 if n_elements(units) eq 0 then units=''
 if n_elements(prj) eq 0 then prj='cyl' 
 if n_elements(sublat) eq 0 then sublat=0.0 
 if n_elements(sublng) eq 0 then sublng=135.0 
 if n_elements(xlab) eq 0 then xlab=''
 if n_elements(ylab) eq 0 then ylab=''
 if n_elements(region) eq 0 then begin 
   if n_elements(minlat) eq 0 then minlat=-90
   if n_elements(maxlat) eq 0 then maxlat=90
   if n_elements(minlon) eq 0 then minlon=0
   if n_elements(maxlon) eq 0 then maxlon=360
 endif 
 if n_elements(nprec) eq 0 then nprec=0
 if n_elements(smooth) eq 0 then smooth=1
 if n_elements(mxt) eq 0 then mxt=0
 if n_elements(myt) eq 0 then myt=0
 if n_elements(misval) eq 0 then misval=-999
 if n_elements(mskval) eq 0 then mskval=-998
 if n_elements(flgval) eq 0 then flgval=-997
 if n_elements(miscol) eq 0 then miscol=0
 if n_elements(mskcol) eq 0 then mskcol=2
 if n_elements(flgcol) eq 0 then flgcol=3
 if n_elements(cth) eq 0 then cth=1
 if n_elements(gth) eq 0 then gth=1
 if n_elements(ath) eq 0 then ath=1
 if n_elements(bcol) eq 0 then bcol=1
 if n_elements(ccol) eq 0 then ccol=0
 if n_elements(gcol) eq 0 then gcol=0
 if n_elements(tcol) eq 0 then tcol=1
 if n_elements(ctyp) eq 0 then ctyp=0 
 if n_elements(gtyp) eq 0 then gtyp=1
 if n_elements(acol) eq 0 then acol=0
 if n_elements(taper) eq 0 then taper=0
 if n_elements(ts) eq 0 then ts=1.0
 if n_elements(ls) eq 0 then ls=1.0
 if n_elements(pos) lt 4 then pos=[0.0,0.0,1.0,1.0]
 if n_elements(cfile) eq 0 then cfile=''
 if n_elements(gfile) eq 0 then gfile='idl.gif'
 if n_elements(pfile) eq 0 then pfile='idl.png'
 if n_elements(psfile) eq 0 then begin 
   if keyword_set(psp) or keyword_set(psl) then psfile='idl.ps'
   if keyword_set(epsp) or keyword_set(epsl) then psfile='idl.eps'
 endif
 if n_elements(xsz) eq 0 then begin
   if keyword_set(psp) or keyword_set(epsp) then xsz=6.5
   if keyword_set(psl) or keyword_set(epsl) then xsz=9.0
 endif
 if n_elements(ysz) eq 0 then begin
   if keyword_set(psp) or keyword_set(epsp) then ysz=9.0
   if keyword_set(psl) or keyword_set(epsl) then ysz=6.5
 endif
 if n_elements(region) eq 4 then begin
   minlat = region(0)
   maxlat = region(1)
   minlon = region(2)
   maxlon = region(3)
   if (minlat gt maxlat) OR (minlon gt maxlon) then begin
     print, 'Error: problem with specified values for region!'
     stop
   endif 
 endif
 if n_elements(aspect) eq 0 then begin
   case prj of 
     'xy':  begin
              aspect=float(maxlon-minlon)/float(maxlat-minlat)
              if (aspect LT 0.1 OR aspect GT 10.0) then $
                aspect=1.0
            end
     'cyl': aspect=float(maxlon-minlon)/float(maxlat-minlat)
     'mer': aspect=float(maxlon-minlon)/float(maxlat-minlat)
     'mol': aspect=float(maxlon-minlon)/float(maxlat-minlat)
     'sat': aspect=1.0
     'polN': begin
               if (maxlon-minlon) gt 180.0 then aspect=1.0 $ 
               else aspect = 1.25
             end
     'polS': begin
               if (maxlon-minlon) gt 180.0 then aspect=1.0 $ 
               else aspect = 1.25
             end
   endcase
 endif
 if n_elements(satpos) eq 0 then satpos=[6.61,0,0] ;GOES elev in globe radii  
 if prj ne 'cyl' AND prj ne 'mer' AND prj ne 'xy' then begin
   left_off=1                      ; No need for rectangular x-y axis 
   bottom_off=1                    ; No need for rectangular x-y axis 
   endif

 
; Open the window or PostScript file or Z buffer for plotting

 ts2=ts
 ls2=ls
 postscript=0
 if keyword_set(psp) OR keyword_set(psl) OR $ 
    keyword_set (epsp) OR keyword_set (epsl) then begin
   postscript=1
   if (!p.font LE 0) then !p.font=0
   set_plot, 'ps'
   load_ctbl, cfile, nc, cmin=cmn, cmax=cmx
   if keyword_set(psp) then begin
     if (!d.unit EQ 0) then device, filename=psfile, /COLOR, BITS_PER_PIXEL=8, $
       /PORTRAIT, xsize=xsz, ysize=ysz, xoffset=1.0, yoffset=1.0, /inches, $
       /Helvetica, font_index=6
   endif 
   if keyword_set(psl) then begin
     if (!d.unit EQ 0) then device, filename=psfile, /COLOR, BITS_PER_PIXEL=8, $
       /LANDSCAPE, xsize=xsz, ysize=ysz, xoffset=1.0, yoffset=10.0, /inches, $
       /Helvetica, font_index=6
   endif
   if keyword_set(epsp) then begin
     if (!d.unit EQ 0) then device, filename=psfile, /COLOR, BITS_PER_PIXEL=8, $
       /PORTRAIT, xsize=xsz, ysize=ysz, xoffset=0.0, yoffset=0.0, /inches, $
       /encapsul, /Helvetica, font_index=6
   endif 
   if keyword_set(epsl) then begin
     if (!d.unit EQ 0) then device, filename=psfile, /COLOR, BITS_PER_PIXEL=8, $
       /LANDSCAPE, xsize=xsz, ysize=ysz, xoffset=0.0, yoffset=9.0, /inches, $
       /encapsul, /Helvetica, font_index=6
   endif
   if (!p.multi(0) EQ 0) then begin
     tvlct, r, g, b, /get
     if (r(bcol) NE 255 OR g(bcol) NE 255 OR b(bcol) NE 255) then begin
       back=bytarr(10,10)+bcol
       TV, back, 0.0, 0.0, XSIZE=1.2, YSIZE=1.2, /norm
     endif
   endif
   xws=!d.x_size
   yws=!d.y_size
   screenflag=1 
 endif 

 if keyword_set(gifout) or keyword_set(pngout) then begin 
   if (!p.multi(0) LE 0) then begin
     set_plot,'z'
     device,set_resolution=[xws,yws], z_buffering=0
   endif else begin
     xws=!d.x_size
     yws=!d.y_size
   endelse
;   if keyword_set(gifout) defsysv, '!GFILE', gfile ; Store output GIF filename
;   if keyword_set(pngout) defsysv, '!PFILE', pfile ; Store output PNG filename
   load_ctbl,cfile,nc,cmin=cmn,cmax=cmx
   if !d.n_colors gt 256 then device, decomposed=0 
   if (!p.multi(0) EQ 0) then erase,bcol
   screenflag=1
   ts2=ts2*(1.0 - (1.0 - float(xws)/600.0)/1.8)
   ls2=ls2*(1.0 - (1.0 - float(xws)/600.0)/1.8)
 endif

 if ( screenflag EQ 0 ) then begin
   set_plot,'x'
   if (!d.window LT 0) then begin
     window,xsize=xws,ysize=yws
   endif else begin
     xws=!d.x_size
     yws=!d.y_size
   endelse
   load_ctbl, cfile, nc, cmin=cmn, cmax=cmx
   if !d.n_colors gt 256 then device, decomposed=0 
   if (!p.multi(0) EQ 0) then erase, bcol
   ts2=ts2*(1.0 - (1.0 - float(xws)/600.0)/1.8)
   ls2=ls2*(1.0 - (1.0 - float(xws)/600.0)/1.8)
 endif 
 
 if keyword_set(multi) then begin
   !p.multi(0)=!p.multi(0)+1
 endif else begin
   !p.multi(0)=0
 endelse
 pmulti=!p.multi(0)
 if keyword_set(order) then !order=1 else !order=0
 
 ; Apply scale and offset to input data; index the missing and masked pixels 

 n=size(a2)
 a3=float(a2)*scale-offset
 if (n[0] NE 2) then begin
   print,'Error, user must input 2d array!'
   stop
 endif 
 dndex=where( a2 NE misval AND a2 NE mskval AND a2 NE flgval, dcnt )
 if (dcnt GT 0) then begin
   if n_elements(vmn) eq 0 then vmn=min(a3(dndex))
   if n_elements(vmx) eq 0 then vmx=max(a3(dndex))
 endif else begin
   if n_elements(vmn) eq 0 then vmn=0
   if n_elements(vmx) eq 0 then vmx=1
 endelse
 mndex = where (a2 EQ misval, mcnt)
 index = where (a2 EQ mskval, lcnt)
 fndex = where (a2 EQ flgval, fcnt)

; Smooth the input data if so requested

 if (keyword_set(smooth)) then begin
   if (smooth GT 1) then begin
     if mcnt GT 0 then a3(mndex) = !VALUES.F_NAN  
     if lcnt GT 0 then a3(index) = !VALUES.F_NAN  
     if fcnt GT 0 then a3(fndex) = !VALUES.F_NAN  
     a3s = smooth ( a3, smooth, /NAN, /edge_truncate ) 
     if mcnt GT 0 then a3s(mndex) = misval
     if lcnt GT 0 then a3s(index) = mskval 
     if fcnt GT 0 then a3s(fndex) = flgval 
     a3 = a3s
   endif
 endif

; Scale the input file to byte values (for use with TV) 

 if (cmx GT nc) then cmx=cmx-(256-nc)

 if (keyword_set(exp)) then begin
   c=expscl( a3, vmn, vmx, cmn, cmx, tcol=tcol, NLEVELS=nct )
   if (vmn EQ -vmx) then begin
     scl=vmx/float(2^(nct/2-1))
     if (n_elements(ctn) EQ 0) then begin
       ctn=strarr(nct+1)
       ctn(nct/2) = fixfmt(0)
       for i=nct/2+1,nct do begin
         x=2^(i-(nct/2+1))
         x=x*scl
         ctn(i) = fixfmt(x)
         ctn(nct-i) = fixfmt(-x)
       endfor
     endif
   endif else begin
     scl=vmx/float(2^nct)
     if (n_elements(ctn) EQ 0) then begin
       ctn=strarr(nct+1)
       for i=0,nct do begin
         x=2^i
         x=x*scl
         ctn(i) = fixfmt(x)
       endfor
     endif
   endelse
 endif else begin
   c=bytscl( a3, MIN=vmn, MAX=vmx, TOP=cmx-cmn ) + cmn
 endelse
 if (mcnt GT 0) then c(mndex)=miscol
 if (lcnt GT 0) then c(index)=mskcol
 if (fcnt GT 0) then c(fndex)=flgcol

; Calculate the position and size of the output image

 imgsize, pos, xws, yws, xs, ys, px, py, xsc, ysc, pxc, pyc, tpos1, tpos2, $
          orient, vmn, vmx, minlon, maxlon, minlat, maxlat, aspect, ts2, ls2, $
          xtl, ytl, ctl, xlab, ylab, title, sttl, units, cbar_off, tick_off, $ 
          left_off, bottom_off, right, top, ytn, vbar, hbar, verbose=verbose
 
 if ( n_elements(imgpos) EQ 0 ) then begin 
   xscl=pos[2]-pos[0]
   yscl=pos[3]-pos[1]
   imgpos=fltarr(8)
   imgpos[0]=(px-pos[0])/xscl
   imgpos[1]=(py-pos[1])/yscl
   imgpos[2]=xs/xscl
   imgpos[3]=ys/yscl
   imgpos[4]=(pxc-pos[0])/xscl
   imgpos[5]=(pyc-pos[1])/yscl
   imgpos[6]=xsc/xscl
   imgpos[7]=ysc/yscl
   if (keyword_set(verbose)) then print,'imgpos = ',imgpos
 endif else begin
   tpos1=tpos1-py-ys
   tpos2=tpos2-py-ys
   xscl=pos[2]-pos[0]
   yscl=pos[3]-pos[1]
   px=imgpos[0]*xscl+pos[0]
   py=imgpos[1]*yscl+pos[1]
   xs=imgpos[2]*xscl
   ys=imgpos[3]*yscl
   pxc=imgpos[4]*xscl+pos[0]
   pyc=imgpos[5]*yscl+pos[1]
   xsc=imgpos[6]*xscl
   ysc=imgpos[7]*yscl
   tpos1=tpos1+py+ys
   tpos2=tpos2+py+ys
 endelse

; Set up map coordinates and projection parameters

 !P.COLOR=acol

 case prj of 
  'cyl':  map_set, 0.0, (minlon+maxlon)/2.0, 0, $ 
                  /cylindrical, limit=[minlat,minlon,maxlat,maxlon], $
                  /noerase, xmargin=[0,0], ymargin=[0,0], $
                  pos=[px,py,px+xs,py+ys], /noborder
  'mer':  begin
            map_set, 0.0, (minlon+maxlon)/2.0, 0, $
            /mercator, limit=[minlat,minlon,maxlat,maxlon], $
            /noerase, xmargin=[0,0], ymargin=[0,0], /clip, $
            pos=[px,py,px+xs,py+ys], /noborder, /label
            ; Projection requires N-S ordering of rows
            if keyword_set(order) then begin
              cnew=reverse(c,2)   
              c=cnew 
              !order=0
            end
          end
  'mol':  begin
            map_set, 0.0, (minlon+maxlon)/2.0, 0, $
            /mollweide, limit=[minlat,minlon,maxlat,maxlon], $
            /noerase, xmargin=[0,0], ymargin=[0,0], $
            pos=[px,py,px+xs,py+ys], /noborder, /label, /horizon 
            ; Projection requires N-S ordering of rows
            if keyword_set(order) then begin
              cnew=reverse(c,2)   
              c=cnew 
              !order=0
            end
          end
  'sat':  begin
            map_set, sublat, sublng, /satellite, $
            limit=[minlat,minlon,maxlat,maxlon], $
            /noerase, xmargin=[0,0], ymargin=[0,0], $
            pos=[px,py,px+xs,py+ys], sat_p=satpos, /noborder, /label
            ; Projection requires N-S ordering of rows
            if keyword_set(order) then begin
              cnew=reverse(c,2)   
              c=cnew 
              !order=0
            end
          end
  'polN': begin
            map_set, 90.0, (minlon+maxlon)/2.0, 0, $
            /stereographic, limit=[0.0, minlon, 90.0, maxlon], $
            /noerase, xmargin=[0,0], ymargin=[0,0], $
            pos=[px,py,px+xs,py+ys], /noborder, /label
            ; Projection requires N-S ordering of rows
            if keyword_set(order) then begin
              cnew=reverse(c,2)   
              c=cnew 
              !order=0
            end
          end
  'polS': begin
            map_set, -90.0, (minlon+maxlon)/2.0, 0, $
            /stereographic, limit=[-90.0, minlon, 0.0, maxlon], $
            /noerase, xmargin=[0,0], ymargin=[0,0], $
            pos=[px,py,px+xs,py+ys], /noborder, /label
            ; Projection requires N-S ordering of rows
            if keyword_set(order) then begin
              cnew=reverse(c,2)   
              c=cnew 
              !order=0
            end
         end
  'xy':  begin
           plot,[minlon,maxlon],[minlat,maxlat],position=[px,py,px+xs,py+ys], $
                xstyle=5,ystyle=5,/nodata,/noerase,/normal
         end
   else: begin
           print, 'Requested projection not recognized' 
           stop         
        end
 endcase

 if (keyword_set(tleft)) then begin 
   xyouts, px, tpos1, string(title,format='("!3",a)'), $
           charsize=ts2*1.5, /normal, alignment=0.0
   xyouts, px, tpos2, string(sttl,format='("!3",a)'), $
           charsize=ts2*0.8, /normal, alignment=0.0
 endif else begin
   xyouts, px+xs/2.0, tpos1, string(title,format='("!3",a)'), $
           charsize=ts2*1.5, /normal, alignment=0.5
   xyouts, px+xs/2.0, tpos2, string(sttl,format='("!3",a)'), $
           charsize=ts2*0.8, /normal, alignment=0.5
 endelse
 !P.CHARSIZE=ls2

; Resize and plot image

 lat1=minlat+0.5*(maxlat-minlat)/float(n(2))
 lat2=maxlat-0.5*(maxlat-minlat)/float(n(2))
 lon1=minlon+0.5*(maxlon-minlon)/float(n(1))
 lon2=maxlon-0.5*(maxlon-minlon)/float(n(1))
 if (keyword_set(verbose)) then print,minlon,maxlon,minlat,maxlat
 if (dcnt GT 0) then begin
   if ( postscript EQ 1 AND (prj EQ 'cyl' OR prj EQ 'xy')) then begin
     TV, c, px, py, XSIZE=xs, YSIZE=ys, /norm
   endif else begin
     if (keyword_set(smooth) AND smooth GT 1) then begin
       if (prj EQ 'xy') then begin
         cc=convert_coord([px,xs],[py,ys],/normal,/to_device)
         px2=cc(0,0)
         py2=cc(1,0)
         xs2=cc(0,1)
         ys2=cc(1,1)
         c2=congrid(c,xs2,ys2)
       endif else begin 
         c2=map_image ( c, px2, py2, xs2, ys2, latmin=lat1, latmax=lat2, $
                        lonmin=lon1, lonmax=lon2, compress=1, /whole_map, $
                        missing=bcol, max_value=cmx )
       endelse
     endif else begin
       if (prj EQ 'xy') then begin
         cc=convert_coord([px,xs],[py,ys],/normal,/to_device)
         px2=cc(0,0)
         py2=cc(1,0)
         xs2=cc(0,1)
         ys2=cc(1,1)
         if (keyword_set(verbose)) then begin
           print,'px  = ',px
           print,'py  = ',py
           print,'xs  = ',xs
           print,'ys  = ',ys
           print,'xs2 = ',xs2
           print,'ys2 = ',ys2
         endif
         c2=congrid(c,xs2,ys2)
       endif else begin 
         c2=map_image ( c, px2, py2, xs2, ys2, latmin=lat1, latmax=lat2, $
                        lonmin=lon1, lonmax=lon2, compress=1, /whole_map, $
                        missing=bcol )
       endelse
     endelse
     TV, c2, px2, py2, XSIZE=xs2, YSIZE=ys2
   endelse
 endif

; Draw the continental outlines and map grid

 if (keyword_set(coast) AND (prj NE 'xy')) then begin
   if keyword_set(usa) then begin
     if keyword_set(islands) then begin
       map_continents, mlinethick=cth, mlinestyle=ctyp, /coasts, /hires, $
                       /usa, color=ccol
     endif else begin
       map_continents, mlinethick=cth, mlinestyle=ctyp, /usa, /continents, $
                       color=ccol
     endelse
   endif else begin
     if keyword_set(islands) then begin
       if (coast EQ 2) then begin
         map_continents, mlinethick=cth, mlinestyle=ctyp, /coasts, /hires, $
                         color=ccol, /fill_continents
       endif else begin
         map_continents, mlinethick=cth, mlinestyle=ctyp, /coasts, /hires, $
                         color=ccol
       endelse
     endif else begin
       if (coast EQ 2) then begin
         map_continents, mlinethick=cth, color=ccol, mlinestyle=ctyp, $
                         /fill_continents
       endif else begin
         map_continents, mlinethick=cth, color=ccol, mlinestyle=ctyp
       endelse
     endelse
   endelse
 endif

; Draw the axis

 case prj of 
  'xy':  xt_axis,nxt,nyt,minlon,maxlon,minlat,maxlat,xtl,ytl,xtn,ytn,xtv,ytv, $
                 xlab,ylab,mxt,myt,left_off,bottom_off,label,right,top, $
                 tlabels,order,thick=ath
  'cyl': map_axis, nxt, nyt, minlon, maxlon, minlat, maxlat, xtl, $
                   ytl, xtn, ytn, xtv, ytv, xlab, ylab, mxt, myt, left_off, $
                   bottom_off, label, right, top, nprec=nprec, thick=ath
  'mer': if keyword_set(label) then begin 
            map_grid, color=gcol, glinethick=gth, glinestyle=gtyp, $
            latdel=(maxlat-minlat)/nyt, londel=(maxlon-minlon)/nxt, /label
         endif else begin
            map_grid, color=gcol, glinethick=gth, glinestyle=gtyp, $
            latdel=(maxlat-minlat)/nyt, londel=(maxlon-minlon)/nxt
         endelse
  'mol': if keyword_set(label) then begin
            map_grid, color=gcol, glinethick=gth, glinestyle=gtyp, $
            latdel=(maxlat-minlat)/nyt, londel=(maxlon-minlon)/nxt, /label
         endif else begin
            map_grid, color=gcol, glinethick=gth, glinestyle=gtyp, $
            latdel=(maxlat-minlat)/nyt, londel=(maxlon-minlon)/nxt
         endelse
  'sat': if keyword_set(label) then begin 
            map_grid, color=gcol, glinethick=gth, glinestyle=gtyp, $
            latdel=(maxlat-minlat)/nyt, londel=(maxlon-minlon)/nxt, /label
         endif else begin
            map_grid, color=gcol, glinethick=gth, glinestyle=gtyp, $
            latdel=(maxlat-minlat)/nyt, londel=(maxlon-minlon)/nxt
         endelse
  'polN': if keyword_set(label) then begin 
             map_grid, color=gcol, glinethick=gth, glinestyle=gtyp, $
             latdel=(maxlat-minlat)/nyt, londel=(maxlon-minlon)/nxt, /label
          endif else begin
             map_grid, color=gcol, glinethick=gth, glinestyle=gtyp, $
             latdel=(maxlat-minlat)/nyt, londel=(maxlon-minlon)/nxt
          endelse
  'polS': if keyword_set(label) then begin 
             map_grid, color=gcol, glinethick=gth, glinestyle=gtyp, $
             latdel=(maxlat-minlat)/nyt, londel=(maxlon-minlon)/nxt, /label
          endif else begin
             map_grid, color=gcol, glinethick=gth, glinestyle=gtyp, $
             latdel=(maxlat-minlat)/nyt, londel=(maxlon-minlon)/nxt
          endelse
 endcase 

; Draw outline around plot

 case prj of 
   'xy':  plot,[minlon,maxlon],[minlat,maxlat],/nodata,/xstyl,/ystyl,xticks=1, $
          xtickn=[' ',' '],yticks=1,ytickn=[' ',' '], $
          position=[px,py,px+xs,py+ys],/noerase,xthick=ath,ythick=ath
   'cyl': plot,[minlon,maxlon],[minlat,maxlat],/nodata,/xstyl,/ystyl,xticks=1, $
          xtickn=[' ',' '],yticks=1,ytickn=[' ',' '], $
          position=[px,py,px+xs,py+ys],/noerase,xthick=ath,ythick=ath
   'mer': plot,[minlon,maxlon],[minlat,maxlat],/nodata,/xstyl,/ystyl,xticks=1, $
          xtickn=[' ',' '],yticks=1,ytickn=[' ',' '], $
          position=[px,py,px+xs,py+ys],/noerase,xthick=ath,ythick=ath
   'mol': map_grid, color=acol,glinestyle=0,glinethick=3,lats=[minlat,maxlat], $
          lons=[minlon,maxlon],xthick=ath,ythick=ath
   'sat': map_horizon, thick=ath, noclip=1
   'polN': map_horizon, thick=ath, noclip=1
   'polS': map_horizon, thick=ath, noclip=1
    else: print,'Plot outline not drawn!'
 endcase

; Draw the colorbar

 !order=0
 if (NOT keyword_set(cbar_off)) then begin
   if (orient EQ 1) then begin
     if (vmn NE vmx) then begin
       colorbar, /vertical, /right, vmin=vmn, vmax=vmx, cmin=cmn, cmax=cmx, $
           ticklen=ctl, position=[pxc,pyc,pxc+xsc,pyc+ysc], title=units, $
           nticks=nct, charsize=ls2, ctickname=ctn, thick=ath, tcol=tcol, $
           taper=taper
     endif
   endif else begin
     if (vmn NE vmx) then begin
       colorbar, /horizontal, /bottom, vmin=vmn, vmax=vmx, cmin=cmn, $
         cmax=cmx, ticklen=ctl, position=[pxc,pyc,pxc+xsc,pyc+ysc], $
         title=units, nticks=nct, charsize=ls2, ctickname=ctn, thick=ath, $
         tcol=tcol, taper=taper
     endif
   endelse
 endif

!p.multi(0)=pmulti

; Close the output PostScript file

 if (keyword_set(psp) OR keyword_set(psl)) AND (!p.multi(0) LE 0) $
 AND NOT (keyword_set(multi)) then begin
   device, /close_file
   set_plot, 'x'
 endif

 if keyword_set(gifout) AND !p.multi(0) LE 0 then begin
  image = tvrd()
  tvlct, /get, r, g, b
  print, gfile
  write_gif, gfile, image, r, g, b
 endif

 if keyword_set(pngout) AND !p.multi(0) LE 0 then begin
  image = tvrd()
  tvlct, /get, r, g, b
  print, pfile
  write_png, pfile, image, r, g, b
 endif

; Set the POSVEC system variable, which makes the page position of the last 
; image drawn available to the next called routine (e.g., overplot.pro)

 pos=[px,py,px+xs,py+ys]
 defsysv, '!POSVEC', exists=i
 if i eq 1 then begin !POSVEC = pos 
 endif else defsysv, '!POSVEC', pos 

; Set the PROJ system variable, which makes the projection of the last image 
; drawn available to the next called routine (e.g., overplot.pro) 

 defsysv, '!PROJ', exists=i
 if i eq 1 then begin !PROJ = prj 
 endif else defsysv, '!PROJ', prj 

; Set the SUBLAT and SUBLNG system variables

 defsysv, '!SUBLAT', exists=i
 if i eq 1 then begin !SUBLAT = sublat 
 endif else defsysv, '!SUBLAT', sublat 

 defsysv, '!SUBLNG', exists=i
 if i eq 1 then begin !SUBLNG = sublng 
 endif else defsysv, '!SUBLNG', sublng 

 defsysv, '!CMIN', exists=i
 if i eq 1 then begin !CMIN = cmn
 endif else defsysv, '!CMIN', cmn

 defsysv, '!CMAX', exists=i
 if i eq 1 then begin !CMAX = cmx
 endif else defsysv, '!CMAX', cmx

; Record the bounding lat/lons in the LIMIT system variable

 defsysv, '!LIMIT', exists=i
 if i eq 1 then begin !LIMIT = [minlat,minlon,maxlat,maxlon] 
 endif else defsysv, '!LIMIT', [minlat,minlon,maxlat,maxlon] 

end