;---------------------------------------------------------------------------
;+
; NAME:
;       PMAP
; PURPOSE:
;       Display or plot a color image of colored dots with axis and a
;       color bar, using one of several projections (default is cylindrical
;       equidistant).
; LAST UPDATE: 11 July 2000 
; CATEGORY:
; CALLING SEQUENCE:
;       pmap      -> Displays help information for pmap.
;       pmap, data, lat, lon
; 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, 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 pmap, data, lat, lon, 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, $
            mis_col=miscol, msk_col=mskcol, 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, fcolor=fcol, 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, 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,' PMAP   -> Displays help information for pmap.'
  print,' PMAP,data,lat,lon -> Displays data 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,'   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,'   PSIZE=psize    Muliplication factor for default pixel size.'
  print,'   ASIZE=asize    Absolute pixel size.'
  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,'   MIS_COL=miscol The color value for missing data. (def=0)'
  print,'   MSK_COL=mskcol The color value for mask data. (def=2)'
  print,'   CCOLOR=ccol    Color for coastal outlines.'
  print,'   FCOLOR=fcol    Color for continental fill.'
  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,x2,y1,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,'   /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), '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), 'ASIZE') gt -1 then asize=xtra.asize
     if strpos(varstr(i-1), 'PSIZE') gt -1 then psize=xtra.psize
     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
   end
 endif

;-------  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 
   ind=where(lat GE  -90.0 AND lat LE  90.0 AND $
             lon GE -180.0 AND lon LE 360.0,cnt)
   if (cnt GT 0) then begin
     if n_elements(minlat) eq 0 then minlat=round(min(lat[ind])-0.5)
     if n_elements(maxlat) eq 0 then maxlat=round(max(lat[ind])+0.5)
     if n_elements(minlon) eq 0 then minlon=round(min(lon[ind])-0.5)
     if n_elements(maxlon) eq 0 then maxlon=round(max(lon[ind])+0.5)
   endif else begin
     print,'No data within bounds, exiting program!'
     stop
   endelse
 endif 
 if n_elements(nprec) eq 0 then nprec=0
 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(miscol) eq 0 then miscol=0
 if n_elements(mskcol) eq 0 then mskcol=2
 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(fcol) eq 0 then fcol=2
 if n_elements(gcol) eq 0 then gcol=0
 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(ts) eq 0 then ts=1.0
 if n_elements(ls) eq 0 then ls=1.0
 if n_elements(psize) eq 0 then psize=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 

 ind=where(lat GE minlat AND lat LE maxlat AND lon GE minlon AND lon LE maxlon,cnt)
 if (cnt GT 0) then begin
   pdata = data[ind]
   plat  = lat[ind]
   plon  = lon[ind]
 endif else begin
   print,'No data within specified region!'
   stop
 endelse

 dndex=where(pdata NE misval AND pdata NE mskval, dcnt )
 if (dcnt GT 0) then begin
   if n_elements(vmn) eq 0 then vmn=min(pdata[dndex])
   if n_elements(vmx) eq 0 then vmx=max(pdata[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 (pdata EQ misval, mcnt)
 index = where (pdata EQ mskval, lcnt)

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

 if (cmx GT nc) then cmx=cmx-(256-nc)
 ncol = cmx - cmn + 1

 if (keyword_set(exp)) then begin
   a=expscl_limits(vmn, vmx, cmn, cmx, NLEVELS=nct )
   lmin=reform(a[0,*])
   lmax=reform(a[1,*])
   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
   lmin=fltarr(ncol)
   lmax=fltarr(ncol)
   for i=0,ncol-1 do begin
     lmin[i]=vmn + float(i)/float(ncol) * (vmx - vmn)
     lmax[i]=vmn + float(i+1)/float(ncol) * (vmx - vmn)
   endfor
 endelse

; 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
          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 
          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
          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
          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
         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*1.0, /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*1.0, /normal, alignment=0.5
 endelse 
 !P.CHARSIZE=ls2

; Determine default symbol size

  if (n_elements(asize) GT 0) then begin
    size = asize
  endif else begin
    npts=n_elements(lat)
    dlat=abs(lat[1:npts-1]-lat[0:npts-2])
    dlon=abs(lon[1:npts-1]-lon[0:npts-2])
    d=sqrt(dlat*dlat+dlon*dlon)
    m=where(d GT 0.0,cnt)
    if (cnt GT 0) then begin
      dsize=min(d[m])
      dx=dsize/(maxlon-minlon)
      dy=dsize/(maxlat-minlat)
    endif else begin
      dx = 0.05
      dy = 0.05
    endelse
    cs=mean([!d.x_ch_size,!d.y_ch_size])
    cx=float(!d.x_size)*(!x.window[1]-!x.window[0])
    cy=float(!d.y_size)*(!y.window[1]-!y.window[0])
    cx=cs/cx
    cy=cs/cy
    size = psize * max([dx/cx,dy/cy])
    if (keyword_set(verbose)) then print,'dx = ',dx
    if (keyword_set(verbose)) then print,'dy = ',dy
    if (keyword_set(verbose)) then print,'cx = ',cx
    if (keyword_set(verbose)) then print,'cy = ',cy
    if (keyword_set(verbose)) then print,'dx/cx = ',dx/cx
    if (keyword_set(verbose)) then print,'dy/cy = ',dy/cy
    if (keyword_set(verbose)) then print,'size = ',size
  endelse
  if (size LT 0.05) then size=0.05
  if (size GT 2.00) then size=2.00
  if (NOT FINITE(size)) then size=0.25
  ;print,'size = ',size
  ;print,'size = ',xwd,yht,csize,dsize,size,format='(a,5(f7.4,x))'
 
; create filled circle for user symbol

  a = findgen(21) * (!PI*2/20)
  usersym,cos(a),sin(a),/fill

; Draw continent map if fill is set

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

; Resize and plot image

 if (keyword_set(verbose)) then print,minlon,maxlon,minlat,maxlat
 m=where(pdata EQ misval,cnt)
 if (cnt GT 0) then oplot, plon[m], plat[m], psym=8, symsize=size,color=miscol
 m=where(pdata EQ mskval,cnt)
 if (cnt GT 0) then oplot, plon[m], plat[m], psym=8, symsize=size,color=mskcol
 if (dcnt GT 0) then begin
   for i=0,ncol-1 do begin
     m=where(pdata GE lmin[i] AND pdata LT lmax[i],cnt)
     ;print,i,lmin[i],lmax[i],cmn+i,cnt
     if (cnt GT 0) then oplot, plon[m], plat[m], psym=8, symsize=size,color=cmn+i
   endfor
   ;m=where(pdata LT lmin[0] AND pdata NE misval AND pdata NE mskval,cnt)
   ;if (cnt GT 0) then oplot, plon[m], plat[m], psym=8, symsize=size,color=cmn
   m=where(pdata GE lmax[ncol-1] AND pdata NE misval AND pdata NE mskval,cnt)
   if (cnt GT 0) then oplot, plon[m], plat[m], psym=8, symsize=size,color=cmx
 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
       map_continents, mlinethick=cth, mlinestyle=ctyp, /coasts, /hires, $
                       color=ccol
     endif else begin
       map_continents, mlinethick=cth, color=ccol, mlinestyle=ctyp
     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
     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
     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