Added the option to provide an output tangent plane for the virtual mosaic. Added an option to provide a reference observation from which to calculate the mosaic's tangent plane

git-svn-id: http://svn.stsci.edu/svn/ssb/stsci_python/stsci_python/trunk/stwcs@11907 fe389314-cf27-0410-b35b-8c050e845b92

1 files changed, 33 insertions, 23 deletions

diff --git a/wcsutil/mosaic.py b/wcsutil/mosaic.py
index 48a8a3a..568b57b 100644
--- a/wcsutil/mosaic.py
+++ b/wcsutil/mosaic.py
@@ -8,7 +8,7 @@ from stwcs.distortion import utils
 from stwcs import updatewcs, wcsutil
 from stwcs.wcsutil import altwcs
 
-def vmosaic(fnames, ext=None, extname=None, undistort=True, wkey='V', wname='VirtualMosaic', plot=False, clobber=False):
+def vmosaic(fnames, outwcs=None, ref_wcs=None, ext=None, extname=None, undistort=True, wkey='V', wname='VirtualMosaic', plot=False, clobber=False):
     """
     Create a virtual mosaic using the WCS of the input images.
     
@@ -16,6 +16,13 @@ def vmosaic(fnames, ext=None, extname=None, undistort=True, wkey='V', wname='Vir
     ----------
     fnames: a string or a list 
               a string or a list of filenames, or a list of wcsutil.HSTWCS objects
+    outwcs: an HSTWCS object
+             if given, represents the output tangent plane
+             if None, the output WCS is calculated from the input observations.
+    ref_wcs: an HSTWCS object
+             if output wcs is not given, this will be used as a reference for the 
+             calculation of the output WCS. If ref_wcs is None and outwcs is None,
+             then the first observation in th einput list is used as reference.
     ext:    an int, a tuple or a list
               an int - represents a FITS extension, e.g. 0 is the primary HDU
               a tuple - uses the notation (extname, extver), e.g. ('sci',1)
@@ -44,27 +51,32 @@ def vmosaic(fnames, ext=None, extname=None, undistort=True, wkey='V', wname='Vir
     Notes
     -----
     The algorithm is:
-    1. The output WCS is calculated based on the input WCSes.
+    1. If output WCS is not given it is calculated from the input WCSes.
        The first image is used as a reference, if no reference is given.
-       This represents the virtual mosaic WCS
+       This represents the virtual mosaic WCS.
     2. For each input observation/chip, an HSTWCS object is created
        and its footprint on the sky is calculated (using only the four corners).
     3. For each input observation the footprint is projected on the output 
        tangent plane and the virtual WCS is recorded in the header.
     """
     wcsobjects = readWCS(fnames, ext, extname)
-    outwcs = utils.output_wcs(wcsobjects, undistort=undistort)
+    if outwcs != None:
+        outwcs = outwcs.deepcopy()
+    else:
+        if ref_wcs != None:
+            outwcs = utils.output_wcs(wcsobjects, ref_wcs=ref_wcs, undistort=undistort)  
+        else:
+            outwcs = utils.output_wcs(wcsobjects, undistort=undistort)
     if plot:
         outc=np.array([[0.,0], [outwcs.naxis1,0], 
                              [outwcs.naxis1, outwcs.naxis2], 
                              [0,outwcs.naxis2], [0,0]])
         plt.plot(outc[:,0], outc[:,1])
-    
     for wobj in wcsobjects:
         outcorners = outwcs.wcs_sky2pix(wobj.calcFootprint(),1)
         if plot:
             plt.plot(outcorners[:,0], outcorners[:,1])
-            objwcs = outwcs.deepcopy()
+        objwcs = outwcs.deepcopy()
         objwcs.wcs.crpix = objwcs.wcs.crpix - (outcorners[0])
         updatehdr(wobj.filename, objwcs,wkey=wkey, wcsname=wname, ext=wobj.extname, clobber=clobber)
     return outwcs
@@ -81,25 +93,26 @@ def updatehdr(fname, wcsobj, wkey, wcsname, ext=1, clobber=False):
             altwcs.deleteWCS(fname, ext=ext, wcskey='V')
     f = pyfits.open(fname, mode='update')
     
-    hwcs = wcs2header(wcsobj, wkey,wcsname).ascardlist()
-    hdr = f[ext].header
-    hdr.ascardlist().extend(hwcs)
+    hwcs = wcs2header(wcsobj)
+    wcsnamekey = 'WCSNAME' + wkey
+    f[ext].header.update(key=wcsnamekey, value=wcsname)
+    for k in hwcs.keys():
+        f[ext].header.update(key=k[:7]+wkey, value=hwcs[k])
+
     f.close()
     
-def wcs2header(wcsobj, wkey='V', wname='VirtualMosaic'):
+def wcs2header(wcsobj):
     
     h = wcsobj.to_header()
-    wcsnamekey = 'WCSNAME' + wkey
-    h.update(key=wcsnamekey, value=wname)
+    
     if wcsobj.wcs.has_cd():
         altwcs.pc2cd(h)
-    for k in h.keys():
-        key = k[:7] + wkey
-        h.update(key=key, value=h[k])
+    h.update('CTYPE1', 'RA---TAN')
+    h.update('CTYPE2', 'DEC--TAN')
     norient = np.rad2deg(np.arctan2(h['CD1_2'],h['CD2_2']))
-    okey = 'ORIENT%s' % wkey
+    #okey = 'ORIENT%s' % wkey
+    okey = 'ORIENT'
     h.update(key=okey, value=norient) 
-    #print h
     return h
 
 def readWCS(input, exts=None, extname=None):
@@ -112,14 +125,11 @@ def readWCS(input, exts=None, extname=None):
                 filelist, output = parseinput.parseinput(input)
             except IOError: raise
     elif isinstance(input, list):
-        if isinstance(input[0], str):
-            # a list of file names
-            try:
-                filelist, output = parseinput.parseinput(input)
-            except IOError: raise
-        elif isinstance(input[0], wcsutil.HSTWCS):
+        if isinstance(input[0], wcsutil.HSTWCS):
             # a list of HSTWCS objects
             return input
+        else:
+            filelist = input[:]
     wcso = []
     fomited = []
     # figure out which FITS extension(s) to use