; PoCaFLIC.PRO ;+ ; NAME: ; POCAFLIC ; PURPOSE: ; Creation of a file containing a series of Polar Camera ; images transformed to global spherical coordinates at a ; fixed height. The user may select either PACE geomagnetic ; coordinates or geographic coordinates. The file is suitable ; for display as a movie using SHOWFLIC.PRO. ; CATEGORY: ; ; CALLING SEQUENCE: ; .RUN POCAFLIC ; INPUTS: ; None. The user is prompted for all needed inputs. ; KEYWORD PARAMETERS: ; None. The user is prompted for all needed inputs. ; OUTPUTS: ; None. ; COMMON BLOCKS: ; HEADBYTES: Contains all header bytes from the images read ; in from the .RAY file, along with exposure start ; UTs, exposure lengths, and flags indicating bad ; images. Shared with RDHBYTES.PRO. ; COORDS: Contains various arrays and parameters used in ; transforming the images to PACE coordinates, and ; displaying the images. Shared with MAP_PGM.PRO or ; MAP_GEO.PRO. ; SIDE EFFECTS: ; Creates a .PGM or .GEO file containing a sequence of 383x383 ; byte images, each one representing a single Polar Camera ; image, mapped to PACE or geographic coordinates, respectively. ; If the user wishes, the images are overlaid with a PACE or ; geographic coordinate grid, and labelled with date, time, ; and altitude information. In the case of geographic ; coordinates, coastlines are also plotted along with the ; coordinate grid. ; RESTRICTIONS: ; ; PROCEDURE: ; ; MODIFICATION HISTORY: ; Written 94/04 by D P Steele. Documented 94/06 by D P Steele. ; Modified 950316 by D P Steele: geographic mapping added as an ; option. ;- COMMON headbytes,oldfname,hbytes,utsec,exptimes,fwtemps,bad_image 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 t0=SYSTIME(1) ; load various string-valued arrays needed @isitdos take2=0B ; Determine graphics device to use. dev='' READ,'Enter desired graphics device ([Z]/X): ',dev dev=STRUPCASE(dev) IF dev NE 'X' THEN dev='Z' !ORDER=0 SET_PLOT,dev IF !D.WINDOW EQ 0 THEN WDELETE,0 ; tell user what imaging modes and filters are available cam_filt,filters ; prompt user and get input cam=0 & filt=cam & day=cam & month=cam PRINT,'Enter the year, month, day, camera, and filter desired' READ,year,month,day,cam,filt PRINT,'Enter the hour, minute, and second to start (default 00:00:00 UT)' get_hms,sthms,shour,sminute,ssecond IF sthms NE '' THEN start=[shour,sminute,ssecond] ELSE start=INTARR(3) PRINT,'Enter the hour, minute, and second to finish (default 23:59:59 UT)' get_hms,fihms,fhour,fminute,fsecond IF fihms NE '' THEN finish=[fhour,fminute,fsecond] ELSE finish=[23,59,59] eheight=[87,100,300,120,100,87,100,120,97,100] PRINT,'Default emission height is '+STRTRIM(eheight(5*cam+filt),2) $ +' km; enter preferred height or to accept this one:' htstr='' READ,htstr IF htstr NE '' THEN eht=ROUND(FLOAT(htstr)) ELSE eht=eheight(5*cam+filt) mtype='' READ,'MLT/MLat dial or geographic coordinates? ([M]/G) ',mtype mtype=STRUPCASE(mtype) IF mtype EQ 'G' THEN BEGIN ctr='' READ,'Center the plot on the geomagnetic pole or Eureka? ([P]/E) ' $ ,ctr pole=(STRUPCASE(ctr) NE 'E') IF NOT pole THEN BEGIN rotstr='' READ,'Rotate the plot to the correct "MLT" orientation? ([Y]/N) ' $ ,rotstr rot=(STRUPCASE(rotstr) NE 'N') ENDIF ELSE rot=1B ENDIF ELSE BEGIN rot=0B pole=0B ENDELSE resp='' READ,'Do you want to label the image? ([Y]/N) ',resp annotate=(STRUPCASE(resp) NE 'N') CASE mtype OF 'M': BEGIN READ,'Do you want to show MLAT/MLT coordinates? ([Y]/N) ',resp grid=(STRUPCASE(resp) NE 'N') END ELSE: BEGIN READ,"Show coastlines ('c') and/or coordinates ('g')? ",resp grid=STRUPCASE(resp) PRINT,'If you want to overplot a satellite footprint, enter' READ,'the standard 2-line element satellite name: ',resp sat=STRUPCASE(resp) END ENDCASE IF (eht LT 300) AND (mtype NE 'G') THEN BEGIN READ,'Do you want to zoom in on the coordinate grid? ([Y]/N)',resp zoom=(STRUPCASE(resp) NE 'N') ENDIF ELSE zoom=0 READ,'Enter the van Rhijn background (R) to subtract: ',vRbg ; Prompt user for fixed saturation brightness PRINT,'If you want to assign the top of the color bar to a fixed number of counts,' PRINT,'enter that number now, otherwise enter anything else for auto-scaling.' tops='' READ,tops flag=isnumber(tops,topnum) ; topnum=0 if tops doesn't represent a number ; ensure a 4-digit year for display; make reasonable assumptions IF year LT 1900 THEN BEGIN IF year LT 50 THEN year=year+2000 IF year GT 50 THEN year=year+1900 ENDIF ; create title for image plot site_names=['','Spyhill Airglow Observatory','Rabbit Lake, SK' $ ,'Eureka AStrO Lab'] PCsite=PoCasite(year,month,day) site=site_names(PCsite) ; Display only those pixels viewing the sky getsig,year,month,day,mask0,mask1 CASE cam OF 0: mask=mask0 1: mask=mask1 ENDCASE ; define path to desired images template=raypath(year,month,day,cam,filt) ; Open the file as a series of small (header-size) units and get ; the header information nimages=rdhbytes(template) ; Load (or create) the data to correct for the reduced transmission. reltrans,oldfname,filtemp,filtrans ; Get the exposure time to be displayed getexpt,exptimes,expt,selexp ; Select images to process todo=getlims(nimages,utsec,start=start,finish=finish) ; Now load the colour table cf=5*cam+filt CASE cf OF 2: ct=3 7: ct=1 8: ct=8 ELSE: ct=0 ENDCASE loadct,ct,/silent ; Open the single big file of processed images and associate it to ; a variable. OPENR,imunit,template,/GET_LUN files=ASSOC(imunit,BYTARR(66048)) ; Open another file to hold the transformed and annotated images. IF mtype EQ 'G' THEN suffix='.geo' ELSE suffix='.pgm' trname=pgmroot+dd $ +STRING(year MOD 100,month,day,cam,filt,suffix $ ,FORMAT='(3I2.2,2I1,A4)') OPENW,trunit,trname,/GET_LUN trans=ASSOC(trunit,BYTARR(383,383)) ; Write information about the mapped image file contents to an ; appropriately-named information file. IF mtype EQ 'G' $ THEN geo_info,template,trname,cam,filt,year,month,day $ ,shour,sminute,ssecond,fhour,fminute,fsecond $ ,eht,vRbg,site,expt,topnum $ ,annotate=annotate,grid=grid,zoom=zoom,sat=sat $ ELSE pgm_info,template,trname,cam,filt,year,month,day $ ,shour,sminute,ssecond,fhour,fminute,fsecond $ ,eht,vRbg,site,expt,topnum $ ,annotate=annotate,zoom=zoom ; Show how long it took to get this far t1=SYSTIME(1) PRINT,t1-t0,FORMAT='(F5.1," s to set up")' ; read in each image, decode the header, rescale from 8 bits, ; transform the image to the desired coordinates, and save it in ; the mapped image file. ip=-1 FOR i=todo(0),todo(1) DO BEGIN whbad=WHERE(bad_image EQ i,nwhbad) IF nwhbad EQ 0 THEN BEGIN KIh=gethd(hbytes(*,i)) ; check exposure time if necessary IF selexp THEN BEGIN thisexp=KIh.misc.exp_scale*KIh.exp.exposure/1000. go_on=(ABS(thisexp-expt) LT 1.0) ENDIF ELSE go_on=1 IF go_on THEN BEGIN file=files(i) im=BYTE(file,512,256,256) im=cleanmod(im,3,2,3,/quiet) ; fix hot pixels ; If file was renormalized to fit 8 bits, reverse this. IF file(290) GT 0 THEN BEGIN IF file(291) EQ 1 THEN im=LONG(im)*(2^file(290)) IF file(291) EQ 0 THEN im=LONG(im)*file(290) ENDIF ; extract filter temperature, determine reduced filter transmission, ; and divide image by reduced transmission to recover actual emission ; brightnesses fwtemp=fwtemps(i) imtrans=interpol(filtrans,filtemp,fwtemp) im=ROUND(im/imtrans(0)) ; Update the image pointer ip=ip+1 ; Set the maximum displayed intensity IF topnum EQ 0 $ THEN lmaxi=pcentile(im*mask,0.999) $ ELSE lmaxi=topnum im(255,98)=lmaxi ; force image to scale properly lmaxbh=pcentile(im*mask,0.001) ; Transform the image to PACE coordinates and label if desired IF mtype EQ 'G' $ THEN map_geo,hbytes(*,i),im*mask,eht,pimg $ ,vRbg=vRbg,fullscale=[lmaxbh,lmaxi] $ ,zoom=zoom,annotate=annotate,grid=grid $ ,rot=rot,pole=pole,sat=sat $ ELSE map_pgm,hbytes(*,i),im*mask,eht,88.664 $ ,-83.444,pimg $ ,vRbg=vRbg,fullscale=[lmaxbh,lmaxi] $ ,zoom=zoom,annotate=annotate ; Save the transformed image to the associated file trans(ip)=pimg EMPTY ; Update the information file with information about the current image pgm_imgs,trname,ip,hbytes(*,i),eur_mlt,plat(asize(1)/2,*) ; Show progress and check for user keypresses PRINT,'.',FORMAT='(A,$)' IF (ip+1) MOD 80 EQ 0 THEN PRINT,FORMAT='()' k=STRUPCASE(GET_KBRD(0)) CASE 1 OF k EQ 'T': BEGIN ; show progress and time of last image PRINT,FORMAT='()' PRINT,todo(0),i,todo(1),FORMAT='(3I4)' PRINT,KIh.misc.tm END k EQ 'X': GOTO,DONE ; force exit RIGHT NOW ELSE: ENDCASE ENDIF ENDIF ENDFOR DONE: CLOSE,imunit & FREE_LUN,imunit CLOSE,trunit & FREE_LUN,trunit PRINT,FORMAT='()' t2=SYSTIME(1) PRINT,t2-t1,ip+1,FORMAT='(F6.1," s to process ",I0," images")' END