PRO map_pgm,headbytes,image,height,eur_plat,eur_plong,pimg $
	   ,new_cam=new_cam,vRbg=vRbg,fullscale=fullscale $
	   ,color=color,zoom=zoom,annotate=annotate

; Keyword FULLSCALE takes 2 values: minimum and maximum intensities
; to be displayed (i.e., limits of colour scale).
; Set keyword COLOR to the color table element you want the coordinate
; grids and annotations drawn in - 0 if you want white.
; Set keyword ZOOM if you want the transformed image to fill the
; plot space optimally; otherwise images mapped to low altitudes
; look small.

COMMON coords,coordht,plat,plong,glat,glong,asize,az,za,col,row $
	     ,vRlayer $
	     ,eur_glat,eur_glong,eur_ht,eur_col,eur_row,eur_mlt $
	     ,mlt_label_offset,plat_contours
COMMON date,year,month,day

; Save the current graphics device and set the device to the Z
; buffer.

olddev=!D.NAME
TVLCT,rx,gx,bx,/GET     ; save the X-windows RGB color table

SET_PLOT,'Z'
DEVICE,Z_BUFFERING=0
loadct,0,/silent
IF N_ELEMENTS(color) EQ 0 THEN color=!D.N_COLORS-1

; If this is the first invocation, read in the coordinate arrays

IF NOT KEYWORD_SET(coordht) THEN coordht=0
IF (height NE coordht) THEN BEGIN
    PRINT,'Reading PGM coordinates for '+STRTRIM(height,2)+' km'
    coordht=FIX(ROUND(height))
    rdpgm,height,plat,plong,glat,glong,/dblsize
    asize=SIZE(plat)

; Zoom in on the coordinates?

    IF KEYWORD_SET(zoom) AND (height LT 300) THEN BEGIN
	oldgrid=FINDGEN(asize(1))
	midgrid=(oldgrid(0)+oldgrid(asize(1)-1))/2.
	newgrid=(oldgrid-midgrid)*(height/300.)+midgrid
	plat=INTERPOLATE(TEMPORARY(plat),newgrid,newgrid,/GRID)
	plong=INTERPOLATE(TEMPORARY(plong),newgrid,newgrid,/GRID)
	glat=INTERPOLATE(TEMPORARY(glat),newgrid,newgrid,/GRID)
	glong=INTERPOLATE(TEMPORARY(glong),newgrid,newgrid,/GRID)
	HELP,plat,plong,glat,glong
	ENDIF

; Define the best plat contours to put on the transformed image plot

    step=nicenumber((90.-MIN(plat(*,asize(1)/2)))/3.)
    plat_contours=90.-[2.,1.]*step

; Define the geographic coordinates of ASTRO Lab

    eur_glat=80.0533
    eur_glong=-86.4161
    eur_ht=0.605

; Find the pixel that best approximates the location of Eureka in
; both geographic and PGM coordinates

    PRINT,'Finding map pixel nearest Eureka'
    find_eur,plat,plong,glat,glong,eur_plat,eur_plong,eur_col,eur_row

; Convert the geographic coordinate arrays to look directions from
; Eureka

    PRINT,'Converting coordinates to look directions'
    geo2az,glat,glong,height,eur_glat,eur_glong,eur_ht,az,za

; If van Rhijn correction is desired, generate the van Rhijn
; background image to be subtracted from all images

    IF KEYWORD_SET(vRbg)  $
	THEN vRlayer=vanRhijn(za,height) $
	ELSE vRlayer=INTARR(256,256)

;Define the ofsets (in degrees) of the MLT labels

    mlt_label_offset=[4.,5.,5.,5.,5.,3.,3.,3.]
    ENDIF

; Display the raw image

; WINDOW,/FREE,XSIZE=asize(1),YSIZE=asize(2),PIXMAP=(N_PARAMS() EQ 6)
; dwin=!D.WINDOW

DEVICE,SET_RESOLUTION=asize(1:2)

chsz=1.5
IF N_ELEMENTS(color) EQ 0 THEN color=!D.N_COLORS-1
lims=pcentile(image,[0.001,0.999])
TVSCL,ROTATE((image>lims(0))<lims(1),2)

; Extract information from the image header

header=gethd(headbytes)
year=header.misc.dt.year
month=header.misc.dt.month
day=header.misc.dt.day
hour=header.misc.tm.hour
minute=header.misc.tm.minute
second=header.misc.tm.second
PCsite=PoCasite(year,month,day)
camera=header.misc.cam
filter=header.misc.fw
exptime=header.misc.exp_scale*header.exp.exposure/1000

IF (N_ELEMENTS(col) EQ 0) OR KEYWORD_SET(new_cam) THEN BEGIN

; Transform the look directions to image coordinates (column,row)

    PRINT,'Converting look directions to column/row'
    az2cr,year,month,day,camera,filter,az,za,col,row $
	 ,/small,/deg,site=PCsite

; Avoid impossible or useless image coordinates

    crej=WHERE((col LE 0) OR (255 LE col),ncrej)
    rrej=WHERE((row LE 0) OR (255 LE row),nrrej)
    IF ncrej GT 0 THEN BEGIN
	col(crej)=0
	row(crej)=0
	ENDIF
    IF nrrej GT 0 THEN BEGIN
	col(rrej)=0
	row(rrej)=0
	ENDIF
    ENDIF

; Sample the mapped image from the observed image

IF N_PARAMS() LT 6 THEN PRINT,'Sampling image'
map=FIX(ROUND((image(col,row)-vRbg*vRlayer(col,row))>0))
gt0=WHERE(map,ngt0)
IF KEYWORD_SET(fullscale)  $
    THEN lims=fullscale $
    ELSE lims=pcentile(map(gt0),[0.001,0.999])
;IF MAX(map) LT lims(1) THEN map(asize(1)/2,asize(2)/2)=lims(1)
map(asize(1)/2,asize(2)/2)=lims(1)	; bright spot at PGM pole

; Get the Eureka MLT values at the beginning of this hour and the next
; hour, and interpolate to find the MLT at the time of the image
IF month GE 10 THEN refyr=year ELSE refyr=year-1
yryr=STRING(refyr MOD 100,(refyr+1) MOD 100,FORMAT='(2I2)')
mlt_file='/jasper/cnsr3_data1/pace/eur_mlt_'+yryr+'.dat'
rdcols,mlt_file,nmlt,yr,mo,dy,hr,mi,se,eurmlt
thishr=WHERE((yr EQ year) AND (mo EQ month) AND $
	     (dy EQ day) AND (hr EQ hour),nthishr)
IF nthishr LT 0 THEN MESSAGE,'No MLT for this hour' ELSE BEGIN
    mlt0=eurmlt(thishr(0))
    mlt1=eurmlt(thishr(0)+1)
    IF mlt1 LT mlt0 THEN mlt1=mlt1+24.
    fhour=(minute*60.+second)/3600.
    eur_mlt=((1.-fhour)*mlt0+fhour*mlt1) MOD 24
    ENDELSE

; Rotate image and maps to give a PGMLT-PGMLat plot

IF N_PARAMS() LT 6 THEN PRINT,'Rotating the mapped image into the MLT-MLat frame'
mlt_rot,map,plat,plong,eur_mlt,eur_col,eur_row,mapcol,maprow

; Display rotated image

c0=0
r0=0
TVSCL,(map>lims(0))<lims(1),c0,r0

; Annotate rotated image

IF KEYWORD_SET(annotate) THEN BEGIN
    PLOT,INDGEN(asize(1)),INDGEN(asize(2)),/DEVICE,/NODATA,/NOERASE $
        ,POSITION=[c0,r0,c0+asize(1),r0+asize(2)],XMARGIN=[0,0],YMARGIN=[0,0] $
        ,XRANGE=[0,asize(1)],YRANGE=[0,asize(2)],XSTYLE=5,YSTYLE=5
    CONTOUR,plat,/OVERPLOT,POSITION=[c0,r0,c0+asize(1),r0+asize(2)],/DEVICE $
           ,XRANGE=[0,asize(1)],YRANGE=[0,asize(2)],XSTYLE=5,YSTYLE=5 $
           ,LEVELS=plat_contours,C_LABELS=[1,1],C_CHARSIZE=chsz,COLOR=color
    rad=asize(1)/2
    FOR i=0,7 DO BEGIN
        ang=-90+i*45
        radii,0,rad,-90+i*45,c0+rad,r0+rad,/DEVICE,COLOR=color
        xoff=170*COS((ang-mlt_label_offset(i))*!DTOR)
        yoff=170*SIN((ang-mlt_label_offset(i))*!DTOR)
        XYOUTS,c0+rad+xoff,r0+rad+yoff,STRTRIM(i*3,2) $
	      ,ALIGNMENT=0.5,/DEVICE,CHARSIZE=chsz,COLOR=color
        IF i EQ 0 THEN XYOUTS,c0+rad-xoff,r0+rad+yoff,'MLT' $
			     ,ALIGNMENT=0.1,CHARSIZE=chsz $
			     ,COLOR=color
        ENDFOR
    XYOUTS,rad,0,STRING(height,FORMAT='(F4.0," km")') $
          ,ALIGNMENT=0.0,/DEVICE,CHARSIZE=chsz,COLOR=color
    PLOTS,c0+mapcol,r0+maprow,PSYM=1,/DEVICE,COLOR=color

    IF KEYWORD_SET(fullscale) THEN BEGIN
        loval=fullscale(0)
        hival=fullscale(1)
        labelimg,headbytes,hival,loval,c0,r0,asize(1),asize(2),chsz $
	        ,color=color,/camfilt,/date
        ENDIF
    ENDIF

EMPTY

; If the mapped image is to be returned, get it from the display
; subsystem

IF (N_PARAMS() EQ 6) THEN pimg=FIX(TVRD()) ELSE k=GET_KBRD(1)
; WDELETE,dwin

DEVICE,/CLOSE

; Restore the graphics device.

SET_PLOT,olddev
TVLCT,rx,gx,bx

RETURN
END