Update to STWCS output_wcs(): This update forces output_wcs() to always return a CD matrix that is orthogonal (cross-terms identical). This is needed to address Ticket #1182.

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

1 files changed, 66 insertions, 2 deletions

diff --git a/lib/stwcs/distortion/utils.py b/lib/stwcs/distortion/utils.py
index 51df391..4228f62 100644
--- a/lib/stwcs/distortion/utils.py
+++ b/lib/stwcs/distortion/utils.py
@@ -1,9 +1,13 @@
 from __future__ import division, print_function # confidence high
 
 import os
+
 import numpy as np
+from numpy import linalg
 from astropy import wcs as pywcs
+
 from stwcs import wcsutil
+from stwcs import updatewcs
 from numpy import sqrt, arctan2
 from stsci.tools import fileutil
 
@@ -36,7 +40,8 @@ def output_wcs(list_of_wcsobj, ref_wcs=None, owcs=None, undistort=True):
         if ref_wcs is None:
             ref_wcs = list_of_wcsobj[0].deepcopy()
         if undistort:
-            outwcs = undistortWCS(ref_wcs)
+            #outwcs = undistortWCS(ref_wcs)
+            outwcs = make_orthogonal_cd(ref_wcs)
         else:
             outwcs = ref_wcs.deepcopy()
         outwcs.wcs.crval = crval
@@ -86,6 +91,66 @@ def computeFootprintCenter(edges):
 
     return crval1,crval2
 
+def make_orthogonal_cd(wcs):
+    """ Create a perfect (square, orthogonal, undistorted) CD matrix from the
+        input WCS.
+    """
+    # get determinant of the CD matrix:
+    cd = wcs.celestial.pixel_scale_matrix
+
+
+    if hasattr(wcs, 'idcv2ref') and wcs.idcv2ref is not None:
+        # Convert the PA_V3 orientation to the orientation at the aperture
+        # This is for the reference chip only - we use this for the
+        # reference tangent plane definition
+        # It has the same orientation as the reference chip
+        pv = updatewcs.makewcs.troll(wcs.pav3,wcs.wcs.crval[1],wcs.idcv2ref,wcs.idcv3ref)
+        # Add the chip rotation angle
+        if wcs.idctheta:
+            pv += wcs.idctheta
+        cs = np.cos(np.deg2rad(pv))
+        sn = np.sin(np.deg2rad(pv))
+        pvmat = np.dot(np.array([[cs,sn],[-sn,cs]]),wcs.parity)
+        rot = np.arctan2(pvmat[0,1],pvmat[1,1])
+        scale = wcs.idcscale/3600.
+
+        det = linalg.det(wcs.parity)
+
+    else:
+
+        det = linalg.det(cd)
+
+        # find pixel scale:
+        if hasattr(wcs, 'idcscale'):
+            scale = (wcs.idcscale) / 3600. # HST pixel scale provided
+        else:
+            scale = np.sqrt(np.abs(det)) # find as sqrt(pixel area)
+
+        # find Y-axis orientation:
+        if hasattr(wcs, 'orientat') and not ignoreHST:
+            rot = np.deg2rad(wcs.orientat) # use HST ORIENTAT
+        else:
+            rot = np.arctan2(wcs.wcs.cd[0,1], wcs.wcs.cd[1,1]) # angle of the Y-axis
+
+    par = -1 if det < 0.0 else 1
+
+    # create a perfectly square, orthogonal WCS
+    sn = np.sin(rot)
+    cs = np.cos(rot)
+    orthogonal_cd = scale * np.array([[par*cs, sn], [-par*sn, cs]])
+
+    lin_wcsobj = pywcs.WCS()
+    lin_wcsobj.wcs.cd = orthogonal_cd
+    lin_wcsobj.wcs.set()
+    lin_wcsobj.orientat = arctan2(lin_wcsobj.wcs.cd[0,1],lin_wcsobj.wcs.cd[1,1]) * 180./np.pi
+    lin_wcsobj.pscale = sqrt(lin_wcsobj.wcs.cd[0,0]**2 + lin_wcsobj.wcs.cd[1,0]**2)*3600.
+    lin_wcsobj.wcs.crval = np.array([0.,0.])
+    lin_wcsobj.wcs.crpix = np.array([0.,0.])
+    lin_wcsobj.wcs.ctype = ['RA---TAN', 'DEC--TAN']
+    lin_wcsobj.wcs.set()
+
+    return lin_wcsobj
+
 def  undistortWCS(wcsobj):
     """
     Creates an undistorted linear WCS by applying the IDCTAB distortion model
@@ -203,4 +268,3 @@ def foundIDCTAB(idctab):
     except KeyError:
         idctab_found = False
     return idctab_found
-