pp_combine implemented

doc/changelog.rst

@@ -3,6 +3,9 @@ Changelog
 
 Major changes to the pipeline since 2016-10-01 (see `Mommert 2017`_) are documented here.
 
+* 2017-10-04: implementation of ``pp_combine``, which enables
+  automated image combination
+
 * 2017-08-29: implementation of Gaia/TGAS as an alternative
   astrometric catalog for shallow widefield observations
 

doc/functions.rst

@@ -393,8 +393,41 @@ Functions that provide additional functionality:
    please refer to the :ref:`manual target identification` walkthrough
    for a recipe on how to use this function.
 
-
+
+.. function:: pp_combine ([-comoving], [-targetname str],
+	      [-manual_rates float, float], [-method str], [-keep_files],
+	      images)
+
+   image combination
+
+   :param comoving: if used, the images will be combined in the moving
+                    frame of a moving target; the target name will be
+                    taken from the ``OBJECT`` header keyword or the
+                    ``targetname`` parameter
+   :param targetname: manual override for the target name if
+                      ``comoving`` parameter is used
+   :param manual_rates: use manual rates instead of queried
+                        ephemerides; in units of arcsec per second in
+                        RA and Dec; RA rate includes factor of cosine
+                        Dec
+   :param method: image combination method: [average, median, clipped]
+                  as provided by `SWARP`_
+   :param keep_files: if used, intermediate files are not deleted
+   :param images: images to run `pp_manident` on
+
+   This function allows the combination of images using different
+   methods. The function makes use of the `SWARP`_ software. By
+   default, images are combined in the rest frame of the background
+   (stars are enhanced, moving objects are partially removed); the
+   ``-comoving`` option enables the combination in the moving frame of
+   one target. In the latter case, images are shifted based on target
+   ephemerides; manual rates can be provided, too. For details on the
+   combination process, please refer to the `SWARP`_ manual. Image
+   files produced by ``pp_combine`` can be used in any other PP
+   function.
+
 
 .. _Source Extractor: http://www.astromatic.net/software/sextractor
 .. _SCAMP: http://www.astromatic.net/software/scamp
 .. _CDS Vizier: http://vizier.u-strasbg.fr/vizier/
+.. _SWARP: http://www.astromatic.net/software/swarp

pp_combine

@@ -0,0 +1 @@
+pp_combine.py

pp_combine.py

@@ -0,0 +1,319 @@
+#!/usr/bin/env python
+
+""" PP_COMBINE - combine frames based on wcs 
+    v1.0: 2017-10-03, michael.mommert@nau.edu
+"""
+from __future__ import print_function
+
+# Photometry Pipeline
+# Copyright (C) 2016  Michael Mommert, michael.mommert@nau.edu
+
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see
+# <http://www.gnu.org/licenses/>.
+
+
+import numpy
+import os
+import sys
+import shutil
+import callhorizons
+import logging
+import subprocess
+import argparse, shlex
+import time
+from astropy.io import fits
+
+### pipeline-specific modules
+import _pp_conf
+import toolbox
+
+# create a portable DEVNULL
+# necessary to prevent subprocess.PIPE and STDOUT from clogging if
+# Source Extractor runs for too long 
+try:
+    from subprocess import DEVNULL # Py3
+except ImportError:
+    import os # Py2
+    DEVNULL = open(os.devnull, 'wb')
+
+# only import if Python3 is used
+if sys.version_info > (3,0):
+    from builtins import str
+
+# setup logging
+logging.basicConfig(filename = _pp_conf.log_filename,
+                    level    = _pp_conf.log_level,
+                    format   = _pp_conf.log_formatline,
+                    datefmt  = _pp_conf.log_datefmt)
+
+
+
+def combine(filenames, comoving, targetname,
+            manual_rate, combine_method, keep_files, display=True,
+            diagnostics=True):
+
+    """
+    image combination wrapper
+    output: diagnostic properties
+    """
+
+    # start logging
+    logging.info('starting image combination with parameters: %s' % \
+                 (', '.join([('%s: %s' % (var, str(val))) for
+                             var, val in list(locals().items())])))
+
+    # check if images have been run through pp_prepare
+    try:
+        midtime_jd = fits.open(filenames[0], verify='silentfix',
+                               ignore_missing_end=True)[0].header['MIDTIMJD']
+    except KeyError:
+        raise KeyError(('%s image header incomplete, have the data run ' +
+                        'through pp_prepare?') % filenames[0])
+        return None
+
+    # adopt first frame as reference frame
+    hdulist = fits.open(filenames[0])
+    header = hdulist[0].header
+    refdate = float(header['MIDTIMJD'])
+    # read out ra and dec from header
+    if obsparam['radec_separator'] == 'XXX':
+        ref_ra_deg = float(header[obsparam['ra']])
+        ref_dec_deg = float(header[obsparam['dec']])    
+        if telescope == 'UKIRTWFCAM':
+            ref_ra_deg = ref_ra_deg/24.*360. - 795/3600.
+            ref_dec_deg -= 795/3600.
+    else:
+        ra_string = header[obsparam['ra']].split(
+            obsparam['radec_separator'])
+        dec_string = header[obsparam['dec']].split(
+            obsparam['radec_separator'])
+        ref_ra_deg = 15.*(float(ra_string[0]) +
+                      old_div(float(ra_string[1]), 60.) +
+                      old_div(float(ra_string[2]), 3600.))
+        ref_dec_deg = (abs(float(dec_string[0])) +
+                   old_div(float(dec_string[1]), 60.) +
+                   old_div(float(dec_string[2]), 3600.))
+        if dec_string[0].find('-') > -1:
+            ref_dec_deg = -1 * dec_deg
+
+    if telescope == 'UKIRTWFCAM':
+        ref_ra_deg = ref_ra_deg/24.*360.
+
+    if obsparam['telescope_keyword'] == "UKIRTWFCAM":
+        ref_ra_deg -= float(header['TRAOFF'])/3600
+        ref_dec_deg -= float(header['TDECOFF'])/3600
+
+    hdulist.close()
+
+    # modify individual frames if comoving == True
+    if comoving:
+        movingfilenames = []
+        for filename in filenames:
+            movingfilename = filename[:filename.find('.fits')]+'_moving.fits'
+            print('shifting %s -> %s' % (filename, movingfilename))
+            logging.info('shifting %s -> %s' % (filename, movingfilename))
+
+            # read out date and pointing information
+            hdulist = fits.open(filename)
+            header = hdulist[0].header
+            date = hdulist[0].header['MIDTIMJD']
+            data = hdulist[0].data
+            hdulist.close()
+
+            # use ephemerides from Horizons if no manual rates are provided
+            if manual_rates is None:
+                # call HORIZONS to get target coordinates
+                eph = callhorizons.query(targetname)
+                eph.set_discreteepochs(date)
+                try:
+                    n = eph.get_ephemerides(str(obsparam['observatory_code']))
+                except ValueError:
+                    print('Target (%s) not an asteroid' % targetname)
+                    logging.warning('Target (%s) not an asteroid' % targetname)
+                    n = None
+
+                if n is None or n == 0:
+                    logging.warning('WARNING: No position from Horizons!'+\
+                                    'Name (%s) correct?' % targetname)
+                    logging.warning('HORIZONS call: %s' % eph.url)
+                    logging.info('proceeding with background sources analysis')
+                    parameters['background_only'] = True
+                else:
+                    logging.info('ephemerides for %s pulled from Horizons' %
+                                 targetname)
+                    target_ra, target_dec = eph[0]['RA'], eph[0]['DEC']
+
+                # get image pointing from header
+                if obsparam['radec_separator'] == 'XXX':
+                    ra_deg = float(header[obsparam['ra']])
+                    dec_deg = float(header[obsparam['dec']])    
+                    if telescope == 'UKIRTWFCAM':
+                        ra_deg = ra_deg/24.*360. - 795/3600.
+                        dec_deg -= 795/3600.
+                else:
+                    ra_string = header[obsparam['ra']].split(
+                        obsparam['radec_separator'])
+                    dec_string = header[obsparam['dec']].split(
+                        obsparam['radec_separator'])
+                    ref_ra_deg = 15.*(float(ra_string[0]) +
+                                      old_div(float(ra_string[1]), 60.) +
+                                      old_div(float(ra_string[2]), 3600.))
+                    ref_dec_deg = (abs(float(dec_string[0])) +
+                                   old_div(float(dec_string[1]), 60.) +
+                                   old_div(float(dec_string[2]), 3600.))
+                    if dec_string[0].find('-') > -1:
+                        ref_dec_deg = -1 * dec_deg
+
+                if filename == filenames[0]:
+                    ref_offset_ra = target_ra - ref_ra_deg
+                    ref_offset_dec = target_dec - ref_dec_deg
+
+                offset_ra = target_ra - ref_ra_deg - ref_offset_ra
+                offset_dec = target_dec - ref_dec_deg - ref_offset_dec
+
+            else:
+                # use manual rates (since they are provided)
+                offset_ra = ((float(header['MIDTIMJD'])-refdate)*86400*
+                             float(manual_rates[0]))/3600
+                offset_dec = ((float(header['MIDTIMJD'])-refdate)*86400*
+                              float(manual_rates[1]))/3600
+
+            logging.info('offsets in RA and Dec: %f, %f arcsec' %
+                         (offset_ra*3600, offset_dec*3600))
+
+            crval1 = float(header['CRVAL1'])
+            crval2 = float(header['CRVAL2']) 
+
+            # write new CRVALi keywords in different file
+            new_hdu = fits.PrimaryHDU(data)
+            new_hdu.header = header
+            new_hdu.header['CRVAL1'] = (crval1-offset_ra,
+                                    'updated in the moving frame of the object')
+            new_hdu.header['CRVAL2'] = (crval2-offset_dec,
+                                    'updated in the moving frame of the object')
+            movingfilenames.append(movingfilename)
+            new_hdu.writeto(movingfilename, overwrite=True,
+                            output_verify='silentfix')
+
+
+
+    if comoving:
+        outfile_name = 'comove.fits'
+        fileline = " ".join(movingfilenames)
+        n_frames = len(movingfilenames)
+    else:
+        outfile_name = 'skycoadd.fits'
+        fileline = " ".join(filenames)        
+        n_frames = len(filenames)
+
+
+    # run swarp on all image catalogs using different catalogs    
+    commandline = (('swarp -combine Y -combine_type %s -delete_tmpfiles '+
+                    'Y -imageout_name %s -interpolate Y -subtract_back N '+
+                    '-weight_type NONE -copy_keywords %s -write_xml N ' +
+                    '-CENTER_TYPE MOST %s') %
+                   ({'median':'MEDIAN', 'average':'AVERAGE',
+                     'clipped':'CLIPPED -CLIP_AMPFRAC 0.2 -CLIP_SIGMA 0.1 '}\
+                    [combine_method],
+                    outfile_name, obsparam['copy_keywords'], fileline))
+
+    logging.info('call SWARP as: %s' % commandline)
+    print('running SWARP to combine {:d} frames...'.format(n_frames), flush=True,
+          end=' ')
+
+    try:
+        swarp = subprocess.Popen(shlex.split(commandline),
+                                 stdout=DEVNULL,
+                                 stderr=DEVNULL,
+                                 close_fds=True)
+        # do not direct stdout to subprocess.PIPE:
+        # for large FITS files, PIPE will clog, stalling 
+        # subprocess.Popen
+    except Exception as e:
+        print('SWARP call:', (e))
+        logging.error('SWARP call:', (e))
+        return None
+
+    swarp.wait()
+    print('done!')
+
+    # remove files that are not needed anymore
+    if not keep_files:
+        if comoving:
+            for filename in movingfilenames:
+                os.remove(filename)
+
+    return n_frames
+
+
+
+
+if __name__ == '__main__':
+
+    # command line arguments
+    parser = argparse.ArgumentParser(description='image combination')
+    parser.add_argument("-comoving", action="store_true",
+                        help='combine in moving target frame')
+    parser.add_argument("-targetname",
+                        help='moving target name')
+    parser.add_argument("-manual_rates", help='manual rates in arcsec/s',
+                        nargs=2)
+    parser.add_argument('-method',
+                        help='combination method',
+                        choices=['average', 'median', 'clipped'],
+                        default='clipped')
+    parser.add_argument("-keep_files", action="store_true",
+                        help='keep intermediate files', default=False)
+    parser.add_argument('images', help='images to process', nargs='+')
+
+    args = parser.parse_args()
+    comoving = args.comoving
+    targetname = args.targetname
+    manual_rates = args.manual_rates
+    combine_method = args.method
+    keep_files = args.keep_files
+    filenames = args.images
+
+    ### read telescope and filter information from fits headers
+    # check that they are the same for all images
+    instruments = []
+    for filename in filenames:
+        hdulist = fits.open(filename, ignore_missing_end=True,
+                            verify='silentfix')
+        header = hdulist[0].header
+        for key in _pp_conf.instrument_keys:
+            if key in header:
+                instruments.append(header[key])
+
+    if len(instruments) == 0:
+        raise KeyError('cannot identify telescope/instrument; please update'
+                       '_pp_conf.instrument_keys accordingly')
+
+    # assign telescope parameters (telescopes.py)
+    telescope = _pp_conf.instrument_identifiers[instruments[0]]
+    obsparam = _pp_conf.telescope_parameters[telescope]
+
+    if manual_rates is not None:
+        comoving = True
+
+    if comoving and targetname is None:
+        targetname = header[obsparam['object']]
+
+    # run image combination wrapper
+    combination = combine(filenames, comoving, targetname,
+                          manual_rates, combine_method, keep_files,
+                          display=True,
+                          diagnostics=True)
+
+