package Subs::UvotImages; ############################################################################## # # DESCRIPTION: Extract transform and plot images # # HISTORY: # HISTORY: $Log: UvotImages.pm,v $ # HISTORY: Revision 1.60 2005/12/22 18:05:09 apsop # HISTORY: Add * for index number in Filename::get method, so that all unfiltered event and raw image files are picked up. # HISTORY: # HISTORY: Revision 1.59 2005/12/16 19:27:44 apsop # HISTORY: Pass attitude to swiftxform when creating grism DET images to get # HISTORY: alternate SKY WCS keys. # HISTORY: # HISTORY: Revision 1.58 2005/12/15 22:37:54 apsop # HISTORY: Pass the expmap flag to uvotimsum to indicate whether or not exposure # HISTORY: maps are being summed. # HISTORY: # HISTORY: Revision 1.57 2005/11/15 13:52:17 apsop # HISTORY: Fix inexplicablly uncaught bug in in mispelling of unlink. # HISTORY: # HISTORY: Revision 1.56 2005/11/09 20:13:56 apsop # HISTORY: Get star catalogue location from sw0.par parameter starcatalog. # HISTORY: # HISTORY: Revision 1.55 2005/11/09 16:30:45 apsop # HISTORY: don;t try and do any image processing if there are no image or event files. # HISTORY: # HISTORY: Revision 1.54 2005/11/08 20:00:18 apsop # HISTORY: <Previous comment bogus>Clean up temporary input list files. # HISTORY: # HISTORY: Revision 1.53 2005/11/08 19:22:28 apsop # HISTORY: Populate the TIMELIST and DATALIST hashes. Used to be an SWCheckInput. # HISTORY: # HISTORY: Revision 1.52 2005/11/03 21:33:01 apsop # HISTORY: Modified a few parameters to work reasonably with Swift Build 16. # HISTORY: # HISTORY: Revision 1.51 2005/11/02 16:34:20 apsop # HISTORY: Use uvotimage to create Level 1 and 2 images. Use uvotskycorr to find # HISTORY: aspect corrections. Only include aspect corrected images in sums. # HISTORY: # HISTORY: Revision 1.50 2005/08/30 14:13:47 apsop # HISTORY: Copy date keywords from sky image prime header to exposure map. # HISTORY: # HISTORY: Revision 1.49 2005/06/01 13:55:27 apsop # HISTORY: Test for presence of event file before trying to use it. # HISTORY: # HISTORY: Revision 1.48 2005/04/29 15:47:24 apsop # HISTORY: Disable calling of uvotmodmap. # HISTORY: # HISTORY: Revision 1.47 2005/04/06 15:43:42 apsop # HISTORY: Change to using CALDB for cal parameters. # HISTORY: # HISTORY: Revision 1.46 2005/03/16 13:31:37 apsop # HISTORY: Add in FILTER keyword to image primary header. # HISTORY: # HISTORY: Revision 1.45 2005/03/15 18:58:39 apsop # HISTORY: Put basic and DATE* keywords into image file primary headers. Fix bug (?) which was not using modmap output for processing. # HISTORY: # HISTORY: Revision 1.44 2005/02/08 18:21:32 apsop # HISTORY: Add new types for uvot filter specific images and exp maps. Replace Filename::corresponding method with file name hack for now. Fix bug which was preventing the production of raw images. # HISTORY: # HISTORY: Revision 1.43 2005/01/21 04:26:43 apsop # HISTORY: Change to use file names stored in cat file rather than contructing them. # HISTORY: # HISTORY: Revision 1.42 2005/01/12 17:29:23 apsop # HISTORY: Change weightfile param in uvotdetect to expfile param. # HISTORY: # HISTORY: Revision 1.41 2005/01/08 04:11:00 apsop # HISTORY: Change test for existence of event list to be more forgiving. # HISTORY: # HISTORY: Revision 1.40 2005/01/07 20:34:21 apsop # HISTORY: No longer necessary to use the CENTER method of swiftxform or a small # HISTORY: value of ncell for uvotmodmap. # HISTORY: # HISTORY: Revision 1.39 2004/12/03 13:38:11 apsop # HISTORY: Turn back on calling uvotmodmap. # HISTORY: # HISTORY: Revision 1.38 2004/10/13 01:38:31 apsop # HISTORY: Disable uvotmodmap for the time being, and fix bug in getting plot file names. # HISTORY: # HISTORY: Revision 1.37 2004/10/12 16:26:36 apsop # HISTORY: Filename changes and a few other tweaks to make module run in complete pipeline. # HISTORY: # HISTORY: Revision 1.36 2004/09/15 22:26:24 apsop # HISTORY: Remove the binning factor from uvot event list file names # HISTORY: # HISTORY: Revision 1.35 2004/09/01 14:42:21 apsop # HISTORY: Disable level 2/3 processing for now. # HISTORY: # HISTORY: Revision 1.34 2004/08/31 20:51:02 apsop # HISTORY: Implemented UVOT level 2 and 3 products. # HISTORY: # HISTORY: Revision 1.33 2004/08/18 17:27:58 apsop # HISTORY: Yet another fix for the timing keywords. # HISTORY: # HISTORY: Revision 1.32 2004/07/19 16:03:24 apsop # HISTORY: Add version number entry to comment field. # HISTORY: # HISTORY: Revision 1.31 2004/07/12 13:46:09 apsop # HISTORY: Fix problems with timing keywords # HISTORY: # HISTORY: Revision 1.30 2004/06/15 16:01:04 apsop # HISTORY: Better trapping of problem uvot exposures. # HISTORY: # HISTORY: Revision 1.29 2004/06/02 16:11:19 apsop # HISTORY: Change name of catalogue extension. # HISTORY: # HISTORY: Revision 1.28 2004/05/04 16:33:01 dah # HISTORY: Bug fixes for case where no initial image file exists. # HISTORY:# HISTORY: Revision 1.27 2004/04/30 16:16:40 dah # HISTORY: Use new FITSfile::keywords() method to speed up image extraction # HISTORY: # HISTORY: Revision 1.26 2004/04/16 20:21:18 dah # HISTORY: Begin using embedded history records # HISTORY: # HISTORY: Revision 1.25 2004/04/16 19:53:18 dah # HISTORY: Begin using embedded history records # HISTORY: # # VERSION: 0.0 # ############################################################################## use Subs::Sub; use Subs::Images; use Util::Xanadu; use Subs::UvotNames; @ISA = ("Subs::Images"); use strict; sub new { my $proto=shift; my $self=$proto->SUPER::new(); $self->{DESCRIPTION}="Extracting, merging and plotting images for UVOT"; return $self; } ################## # METHODS: ################## sub body { my $self=shift; $self->create_images(); $self->image_raw_to_det(); $self->aspect_correct_sky_images(); $self->create_exposure_maps(); # create summed images/exposure maps $self->sum_images(); # detect sources and determine magnitudes $self->detect_sources(); $self->plot_images(); } # end of body method ############################################################################# # Extract images from filtered event data and raw images files ############################################################################# sub create_images { my $self = shift; my $log =$self->log(); my $filename=$self->filename(); my $procpar =$self->procpar(); my $jobpar =$self->jobpar(); $log->entry('Creating images from UVOT event and image mode data'); my $attfile = $filename->get('attitude', 's'); if (not $attfile) { $log->error(2, 'not attitude file available'); return; } my @eventFiles = $filename->get('unfiltered', 'uvot', '*', '*'); my @imageFiles = $filename->get('rawimage', 'uvot', '*', '*'); unless( @eventFiles || @imageFiles ){ $log->error(1, 'No image or event files to process.'); return; } my $infile = 'uvotimage.files'; my $fh = FileHandle->new($infile, 'w'); if (not $fh) { $log->error(2, "unable to create $infile [$!]"); return; } foreach my $name (@eventFiles, @imageFiles) { $fh->print($name . "\n"); } $fh->close; my $catfile = $filename->get('hk', 'uvot', 'ct', '*'); if (not -f $catfile) { $log->entry('No UVOT exposure catalogue to update'); } else { $log->entry('UVOT exposure catalogue update not implemented'); } my $prefix = 'Qz3'; # a unique string my $tool = Util::HEAdas->new('uvotimage')->is_script(1); $tool->params({ infile => '@' . $infile, prefix => $prefix, attfile => $attfile, teldeffile => 'CALDB', alignfile => 'CALDB', ra => $jobpar->read('ra'), dec => $jobpar->read('dec'), roll => $jobpar->read('roll'), # badpix => 'no', # flatfield => 'no', # mod8corr => 'no', # catfile => $catfile, chatter => 5, }) ->run; unlink $infile; =pod # edit each image to include keywords my $dateobs = $jobpar->read('obsdate') . 'T' . $jobpar->read('obstime'); my $dateend = $jobpar->read('enddate') . 'T' . $jobpar->read('endtime'); my %keys = ( 'DATE-OBS' => $dateobs, 'DATE-END' => $dateend, RA_PNT => $jobpar->read('burst_ra'), DEC_PNT => $jobpar->read('burst_dec'), PA_PNT => $jobpar->read('burst_roll'), ); my $keywords = 'uvotimage.key'; $fh = FileHandle->new($keywords, 'w'); while (my ($k, $v) = each(%keys)) { $fh->print("$k = $v\n"); } $fh->close; my $edit = Util::HEAdas->new('fthedit') ->params({ keywords => '@' . $keywords, }); =cut # move each file ${prefix}xxx to sw<obsid>xxx my $swobsid = $filename->sequence_specific; foreach my $name (glob($prefix . '*')) { my $xxx = substr($name, length($prefix)); my $fixed = $swobsid . $xxx; if (-f $fixed) { unlink($fixed); } rename($name, $fixed); # $edit->params({ infile => $fixed })->run; } # unlink($keywords); } # end create images method ############################################################################### # convert raw coordinate images to detector coordinates ############################################################################### sub image_raw_to_det { my $self=shift; my $log =$self->log(); my $filename=$self->filename(); my $procpar =$self->procpar(); my $jobpar =$self->jobpar(); $log->entry("Converting raw grism images to detector coordinates"); my $attfile = $filename->get('attitude', 's'); if (not $attfile) { $attfile = 'NONE'; } ######################################### # create the swiftxform tool ######################################### my $swiftxform = Util::HEAdas->new('swiftxform'); $swiftxform->params({ teldeffile => 'CALDB', to => 'DET', attfile => $attfile, ra => $jobpar->read('ra'), dec => $jobpar->read('dec'), roll => $jobpar->read('roll'), aberration => 'no', seed => $procpar->read('seed'), chatter => 4, }) ->is_script(1); ###################### # loop over _GRISM_ files # since only they are converted to DET coordinates ###################### foreach my $rawFile ($filename->get('rawimage', 'uvot', 'g*', '*') ) { ######################################### # determine the DET coordinate filename ######################################### my $detFile = $filename->corresponding("rawimage", "detimage", $rawFile); $log->entry("converting $rawFile to $detFile"); ###################### # do the conversion ###################### $swiftxform->params({ infile => $rawFile, outfile => $detFile, }) ->run(); } # end of loop over files } # end of image_raw_to_det method sub aspect_correct_sky_images { my $self = shift; my $log =$self->log(); my $filename=$self->filename(); my $procpar =$self->procpar(); my $jobpar =$self->jobpar(); $log->entry('Performing aspect correction on UVOT images'); my $attfile = $filename->get('attitude', 's'); if (not $attfile) { $log->error(2, 'no attitude file available'); return; } my @skyFiles = $filename->get('filterimg', 'uvot', '*'); if (not @skyFiles) { $log->entry("no sky images to correct"); return; } my $infile = 'uvotskycorr.files'; my $fh = FileHandle->new($infile, 'w'); if (not $fh) { $log->error(2, "unable to create $infile [$!]"); return; } foreach my $name (@skyFiles) { $fh->print($name . "\n"); } $fh->close; my $corrfile = 'ASPCORR.txt'; my $catspec = $procpar->read('starcatalog'); unlink($corrfile); $log->entry('finding corrections'); my $find = Util::HEAdas->new('uvotskycorr')->is_script(1); $find->params({ skyfile => '@' . $infile, what => 'ID', outfile => $corrfile, corrfile => 'NONE', attfile => $attfile, partition => $catspec, starid => 'n.reference=50 n.observation=30', chatter => 5, }) ->run; return if not -f $corrfile; $log->entry('applying corrections'); my $apply = Util::HEAdas->new('uvotskycorr')->is_script(1); $apply->params({ skyfile => '@' . $infile, what => 'SKY', outfile => 'NONE', corrfile => $corrfile, attfile => $attfile, partition => $catspec, chatter => 5, }) ->run; unlink $infile; } # end of aspect_correct_sky_images method ############################################################################### # generate exposure maps ############################################################################### sub create_exposure_maps { my $self=shift; my $log = $self->log(); my $filename= $self->filename(); my $procpar = $self->procpar(); my $jobpar = $self->jobpar(); $log->entry("Creating exposure maps"); ###################################################### # make sure there is an attitude file ###################################################### my $attitude = $filename->get('attitude', 's'); unless(-f $attitude) { $log->entry("No attitude data available - can't make exposure maps"); return; } ######################################### # create the swiftxform tool ######################################### my $uexpmap = Util::HEAdas->new('uvotexpmap'); $uexpmap->params({ attfile => $attitude, teldeffile => 'CALDB', method => 'MEANFOV', # TODO: is SHIFT/ADD in effect? attdelta => 100, # TODO: make a parameter aberration => 'no', chatter => 4, }); ###################### # loop over files ###################### foreach my $skyFile ($filename->get('filterimg', 'uvot', '*', '*') ) { # exclude grism files next if $filename->is_grism($skyFile, 'skyimage'); ######################################### # determine the exposure map filename ######################################### my $expFile = $skyFile; $expFile =~ s/_sk\./_ex\./; $log->entry("building exposure map $expFile for $skyFile"); my $badFile = $skyFile; $badFile =~ s/_sk\./_bp\./; if (not -f $badFile) { $log->entry("no bad pixel file for $skyFile"); $badFile = 'NONE'; } ####################### # make the exposure map ####################### $uexpmap->params({ infile => $skyFile, badpixfile => $badFile, outfile => $expFile, }) ->is_script(1) ->run(); } # end of loop over files } # end of create_exposure_maps method sub sum_images { my $self=shift; my $log = $self->log(); my $filename= $self->filename(); my $procpar = $self->procpar(); my $jobpar = $self->jobpar(); $log->entry("Creating summed images"); ################################ # Create tool for summing images ################################ my $usumexp = Util::HEAdas->new('uvotimsum'); my $tmpFile = 'sum.tmp'; ######################################## # clean up any pre-existing summed files ######################################## foreach my $type (qw(skyimage expimage)) { foreach my $imageFile ($filename->get($type, 'uvot', '', '*')) { $log->entry("removing obsolete summed image $imageFile"); unlink($imageFile); } } ###################### # loop over files ###################### my %filter_to_sum = ('filterimg' => 'skyimage', 'filterexp' => 'expimage'); foreach my $type (keys %filter_to_sum) { foreach my $imageFile ($filename->get($type, 'uvot', '*', '*')) { my (undef, $filter_code) = $filename->parse($imageFile, $type); $log->entry("Filter code for $imageFile is $filter_code."); next if not $filter_code; ######################################### # determine the summed image filename # by removing the filter code from sw\d{9}uff(sk|ex).img ######################################### my $sumFile = $filename->get($filter_to_sum{$type}, 'uvot', '', 0); $log->entry("summing $imageFile to $sumFile"); ###################################################### # determine which extensions to include in the sum # (i.e., which extensions have been aspect corrected) ###################################################### my $nhdu = 0; my $fits = Util::SimpleFITS->readonly($imageFile) ->nhdu($nhdu); my $status = 0; if ($status = $fits->status) { $log->error(1, "unable to open $imageFile [$status]"); next; } my %include; my %exclude; for (my $i = 2; $i <= $nhdu; ++$i) { my $extname = ''; my $aspcorr = ''; my $tmp = 0; $fits->move($i) ->readkey(EXTNAME => $extname) ->readkey(ASPCORR => $aspcorr) ->status($tmp) ->setstatus(0) ; if ($tmp) { $log->entry("unable to read EXTNAME/ASPCORR from HDU $i [$tmp]"); } elsif ($aspcorr =~ /DIRECT/) { $include{$extname} = 1; } if (not $include{$extname}) { $exclude{$extname} = 1; } } if ($status = $fits->close->status) { $log->error(1, "unable to close $imageFile [$status]"); } # only create a summed image if there are aspect corrected images if (not scalar(keys(%include))) { $log->entry("no aspect corrected images in $imageFile to sum"); next; } my $exclude = join(',', keys(%exclude)); if ($exclude) { $log->entry("excluding extensions $exclude from sum"); } my $filter = $filename->filter_name_for_code($filter_code); ###################### # do the summation ###################### $usumexp->params({ infile => $imageFile, outfile => $tmpFile, exclude => ($exclude || 'NONE'), expmap => ($type =~ /exp/ ? 'yes' : 'no'), }) ->is_script(1) ->run(); my $keys = "[col " . " #EXTNAME='$filter';" . " #FILTER='$filter';" . ']'; if (not -e $sumFile) { Util::HEAdas->new('ftcopy') ->params({ infile => $tmpFile . "[col #EXTNAME='$filter']", outfile => $sumFile, copyall => 'yes', }) ->run(); } else { # append this result to the sum file Util::HEAdas->new('ftappend') ->params({ infile => "$tmpFile+1[col #EXTNAME='$filter']", outfile => $sumFile, }) ->run(); } unlink($tmpFile); } # end of loop over files } # end of loop over types } # end of process_sky method sub detect_sources { my $self=shift; my $log = $self->log(); my $filename= $self->filename(); my $procpar = $self->procpar(); my $jobpar = $self->jobpar(); $log->entry("Detecting sources"); ######################################## # clean up any pre-existing source lists ######################################## foreach my $catFile ($filename->get('srclist', 'uvot', '', '*')) { unlink($catFile); } ############################################################### # Create tools for detecting sources and determining magnitudes ############################################################### my $udetect = Util::HEAdas->new('uvotdetect'); my $umag = Util::HEAdas->new('uvotmag'); $umag->params({ zerofile => 'CALDB', coinfile => 'CALDB', filter => 'default', }); my $tmpFile = 'sources.tmp'; my $catFile = undef; ############################# # loop over summed sky images ############################# foreach my $skyFile ($filename->get('skyimage', 'uvot', '', '*')) { $log->entry("detect: skyFile '$skyFile'"); my $expFile = $filename->corresponding('skyimage', 'expimage', $skyFile); if (not $expFile or not -f $expFile) { $log->entry("No exposure data for $skyFile - assuming constant exposure"); $expFile = 'NONE'; } if (not $catFile) { $catFile = $filename->corresponding('skyimage', 'srclist', $skyFile); } ######################################### # run uvotdetect on each extension ######################################### my $fitsFile = Util::FITSfile->new($skyFile, 0); my $nhdus = $fitsFile->nhdus(); for (my $hdu0 = 1; $hdu0 < $nhdus; ++$hdu0) { $fitsFile->ext($hdu0); my $filter = $fitsFile->keyword('EXTNAME'); my $weightFile = $expFile; if ($expFile ne 'NONE') { $weightFile .= "[$filter]"; } if (-e $tmpFile) { unlink($tmpFile); } $udetect->params({ infile => "$skyFile+$hdu0", outfile => $tmpFile, expfile => $weightFile, zerobkg => -1, threshold => 5, # TODO: make parameter }) ->is_script(1) ->run(); # run uvotmag if (-e $tmpFile) { $umag->params({ infile => $tmpFile, }) ->is_script(1) ->run(); } if (not -e $catFile) { Util::HEAdas->new('ftcopy') ->params({ infile => $tmpFile . "[col #EXTNAME='$filter']", outfile => $catFile, copyall => 'yes', }) ->run(); } else { # append this result to the sum file Util::HEAdas->new('ftappend') ->params({ infile => "$tmpFile+1[col #EXTNAME='$filter']", outfile => $catFile, }) ->run(); } } # end of loop over extensions } # end of loop over files unlink($tmpFile) if -e $tmpFile; } # end of detect_sources method ################################################################################ # ################################################################################ sub plot_images { my $self = shift; my $log =$self->log(); my $filename=$self->filename(); my $procpar =$self->procpar(); my $jobpar =$self->jobpar(); my $max_image_dimen=1024; $log->entry("Plotting all images"); my $ximage = Util::Xanadu->new("ximage"); my $type; foreach $type ("detimage", "skyimage") { my $plot_type = $type; $plot_type =~ s/image/plot/; ###################################### # loop over summed images of this type ###################################### my $file; foreach my $file ($filename->get($type, 'uvot', '', '*')) { my $fitsFile = Util::FITSfile->new($file, 0); my $nhdus = $fitsFile->nhdus(); my $grid=''; if($type =~ /^sky/ ) { $grid = 'grid' } ############################################## # for each extension (filter) create a graphic ############################################## for(my $hdu=1; $hdu<$nhdus; $hdu++) { my @commands; push @commands, ('cpd /vgif', 'cey 2000'); $fitsFile->ext($hdu); my $filter = $fitsFile->keyword('EXTNAME'); my $filter_code = $filename->filter_code_for_name($filter); my $plotFile = $filename->get($plot_type, 'uvot', $filter_code, 0); ######################################## # determine if we need to bin the image ######################################## my $dimen = $fitsFile->keyword('NAXIS1'); my $bin=1; for($bin=1; $max_image_dimen*$bin<$dimen; $bin++) {} $log->entry("Plotting ${file}\[$hdu\] in $plotFile binned by $bin"); push @commands, ("read/fits/rebin=$bin \{${file}\[$hdu\]\}", 'smooth', 'disp', $grid, 'scale'); push @commands, 'exit'; $ximage->script(@commands)->run(); rename 'pgplot.gif', $plotFile; } } } } # end of plot_images method 1;