Merge pull request #9 from nden/refactor-and-tests

restructure and add stwcs tests

1 files changed, 988 insertions, 0 deletions

diff --git a/stwcs/wcsutil/hstwcs.py b/stwcs/wcsutil/hstwcs.py
new file mode 100644
index 0000000..bfebcfc
--- /dev/null
+++ b/stwcs/wcsutil/hstwcs.py
@@ -0,0 +1,988 @@
+from __future__ import absolute_import, division, print_function # confidence high
+
+import os
+from astropy.wcs import WCS
+from astropy.io import fits
+from stwcs.distortion import models, coeff_converter
+import numpy as np
+from stsci.tools import fileutil
+
+from . import altwcs
+from . import getinput
+from . import mappings
+from . import instruments
+from .mappings import inst_mappings, ins_spec_kw
+from .mappings import basic_wcs
+
+__docformat__ = 'restructuredtext'
+
+#
+#### Utility functions copied from 'updatewcs.utils' to avoid circular imports
+#
+def extract_rootname(kwvalue,suffix=""):
+    """ Returns the rootname from a full reference filename
+
+        If a non-valid value (any of ['','N/A','NONE','INDEF',None]) is input,
+            simply return a string value of 'NONE'
+
+        This function will also replace any 'suffix' specified with a blank.
+    """
+    # check to see whether a valid kwvalue has been provided as input
+    if kwvalue.strip() in ['','N/A','NONE','INDEF',None]:
+        return 'NONE' # no valid value, so return 'NONE'
+
+    # for a valid kwvalue, parse out the rootname
+    # strip off any environment variable from input filename, if any are given
+    if '$' in kwvalue:
+        fullval = kwvalue[kwvalue.find('$')+1:]
+    else:
+        fullval = kwvalue
+    # Extract filename without path from kwvalue
+    fname = os.path.basename(fullval).strip()
+
+    # Now, rip out just the rootname from the full filename
+    rootname = fileutil.buildNewRootname(fname)
+
+    # Now, remove any known suffix from rootname
+    rootname = rootname.replace(suffix,'')
+    return rootname.strip()
+
+def build_default_wcsname(idctab):
+
+    idcname = extract_rootname(idctab,suffix='_idc')
+    wcsname = 'IDC_' + idcname
+    return wcsname
+
+
+class NoConvergence(Exception):
+    """
+    An error class used to report non-convergence and/or divergence of
+    numerical methods. It is used to report errors in the iterative solution
+    used by the :py:meth:`~stwcs.hstwcs.HSTWCS.all_world2pix`\ .
+
+    Attributes
+    ----------
+
+    best_solution : numpy.array
+        Best solution achieved by the method.
+
+    accuracy : float
+        Accuracy of the :py:attr:`best_solution`\ .
+
+    niter : int
+        Number of iterations performed by the numerical method to compute
+        :py:attr:`best_solution`\ .
+
+    divergent : None, numpy.array
+        Indices of the points in :py:attr:`best_solution` array for which the
+        solution appears to be divergent. If the solution does not diverge,
+        `divergent` will be set to `None`.
+
+    failed2converge : None, numpy.array
+        Indices of the points in :py:attr:`best_solution` array for which the
+        solution failed to converge within the specified maximum number
+        of iterations. If there are no non-converging poits (i.e., if
+        the required accuracy has been achieved for all points) then
+        `failed2converge` will be set to `None`.
+
+    """
+    def __init__(self, *args, **kwargs):
+        super(NoConvergence, self).__init__(*args)
+
+        self.best_solution  = kwargs.pop('best_solution', None)
+        self.accuracy       = kwargs.pop('accuracy', None)
+        self.niter          = kwargs.pop('niter', None)
+        self.divergent      = kwargs.pop('divergent', None)
+        self.failed2converge= kwargs.pop('failed2converge', None)
+
+
+#
+#### HSTWCS Class definition
+#
+class HSTWCS(WCS):
+
+    def __init__(self, fobj=None, ext=None, minerr=0.0, wcskey=" "):
+        """
+        Create a WCS object based on the instrument.
+
+        In addition to basic WCS keywords this class provides
+        instrument specific information needed in distortion computation.
+
+        Parameters
+        ----------
+        fobj : str or `astropy.io.fits.HDUList` object or None
+            file name, e.g j9irw4b1q_flt.fits
+            fully qualified filename[EXTNAME,EXTNUM], e.g. j9irw4b1q_flt.fits[sci,1]
+            `astropy.io.fits` file object, e.g fits.open('j9irw4b1q_flt.fits'), in which case the
+            user is responsible for closing the file object.
+        ext : int, tuple or None
+            extension number
+            if ext is tuple, it must be ("EXTNAME", EXTNUM), e.g. ("SCI", 2)
+            if ext is None, it is assumed the data is in the primary hdu
+        minerr : float
+            minimum value a distortion correction must have in order to be applied.
+            If CPERRja, CQERRja are smaller than minerr, the corersponding
+            distortion is not applied.
+        wcskey : str
+            A one character A-Z or " " used to retrieve and define an
+            alternate WCS description.
+        """
+
+        self.inst_kw = ins_spec_kw
+        self.minerr = minerr
+        self.wcskey = wcskey
+
+        if fobj is not None:
+            filename, hdr0, ehdr, phdu = getinput.parseSingleInput(f=fobj,
+                                                                   ext=ext)
+            self.filename = filename
+            instrument_name = hdr0.get('INSTRUME', 'DEFAULT')
+            if instrument_name == 'DEFAULT' or instrument_name not in list(inst_mappings.keys()):
+                #['IRAF/ARTDATA','',' ','N/A']:
+                self.instrument = 'DEFAULT'
+            else:
+                self.instrument = instrument_name
+            # Set the correct reference frame
+            refframe = determine_refframe(hdr0)
+            ehdr['RADESYS'] = refframe
+
+            WCS.__init__(self, ehdr, fobj=phdu, minerr=self.minerr,
+                         key=self.wcskey)
+            if self.instrument == 'DEFAULT':
+                self.pc2cd()
+            # If input was a `astropy.io.fits.HDUList` object, it's the user's
+            # responsibility to close it, otherwise, it's closed here.
+            if not isinstance(fobj, fits.HDUList):
+                phdu.close()
+            self.setInstrSpecKw(hdr0, ehdr)
+            self.readIDCCoeffs(ehdr)
+            extname = ehdr.get('EXTNAME', '')
+            extnum = ehdr.get('EXTVER', None)
+            self.extname = (extname, extnum)
+        else:
+            # create a default HSTWCS object
+            self.instrument = 'DEFAULT'
+            WCS.__init__(self, minerr=self.minerr, key=self.wcskey)
+            self.pc2cd()
+            self.setInstrSpecKw()
+        self.setPscale()
+        self.setOrient()
+
+    @property
+    def naxis1(self):
+        return self._naxis1
+
+    @naxis1.setter
+    def naxis1(self, value):
+        self._naxis1 = value
+
+    @property
+    def naxis2(self):
+        return self._naxis2
+
+    @naxis2.setter
+    def naxis2(self, value):
+        self._naxis2 = value
+
+    def readIDCCoeffs(self, header):
+        """
+        Reads in first order IDCTAB coefficients if present in the header
+        """
+        coeffs = ['ocx10', 'ocx11', 'ocy10', 'ocy11', 'idcscale',
+                    'idcv2ref','idcv3ref', 'idctheta']
+        for c in coeffs:
+            self.__setattr__(c, header.get(c, None))
+
+    def setInstrSpecKw(self, prim_hdr=None, ext_hdr=None):
+        """
+        Populate the instrument specific attributes:
+
+        These can be in different headers but each instrument class has knowledge
+        of where to look for them.
+
+        Parameters
+        ----------
+        prim_hdr : `astropy.io.fits.Header`
+            primary header
+        ext_hdr : `astropy.io.fits.Header`
+            extension header
+
+        """
+        if self.instrument in list(inst_mappings.keys()):
+            inst_kl = inst_mappings[self.instrument]
+            inst_kl = instruments.__dict__[inst_kl]
+            insobj = inst_kl(prim_hdr, ext_hdr)
+
+            for key in self.inst_kw:
+                try:
+                    self.__setattr__(key, insobj.__getattribute__(key))
+                except AttributeError:
+                    # Some of the instrument's attributes are recorded in the primary header and
+                    # were already set, (e.g. 'DETECTOR'), the code below is a check for that case.
+                    if not self.__getattribute__(key):
+                        raise
+                    else:
+                        pass
+
+        else:
+            raise KeyError("Unsupported instrument - %s" %self.instrument)
+
+    def setPscale(self):
+        """
+        Calculates the plate scale from the CD matrix
+        """
+        try:
+            cd11 = self.wcs.cd[0][0]
+            cd21 = self.wcs.cd[1][0]
+            self.pscale = np.sqrt(np.power(cd11,2)+np.power(cd21,2)) * 3600.
+        except AttributeError:
+            if self.wcs.has_cd():
+                print("This file has a PC matrix. You may want to convert it \n \
+                to a CD matrix, if reasonable, by running pc2.cd() method.\n \
+                The plate scale can be set then by calling setPscale() method.\n")
+            self.pscale = None
+
+    def setOrient(self):
+        """
+        Computes ORIENTAT from the CD matrix
+        """
+        try:
+            cd12 = self.wcs.cd[0][1]
+            cd22 = self.wcs.cd[1][1]
+            self.orientat = np.rad2deg(np.arctan2(cd12,cd22))
+        except AttributeError:
+            if self.wcs.has_cd():
+                print("This file has a PC matrix. You may want to convert it \n \
+                to a CD matrix, if reasonable, by running pc2.cd() method.\n \
+                The orientation can be set then by calling setOrient() method.\n")
+            self.pscale = None
+
+    def updatePscale(self, scale):
+        """
+        Updates the CD matrix with a new plate scale
+        """
+        self.wcs.cd = self.wcs.cd/self.pscale*scale
+        self.setPscale()
+
+    def readModel(self, update=False, header=None):
+        """
+        Reads distortion model from IDCTAB.
+
+        If IDCTAB is not found ('N/A', "", or not found on disk), then
+        if SIP coefficients and first order IDCTAB coefficients are present
+        in the header, restore the idcmodel from the header.
+        If not - assign None to self.idcmodel.
+
+        Parameters
+        ----------
+        header : `astropy.io.fits.Header`
+            fits extension header
+        update : bool (False)
+            if True - record the following IDCTAB quantities as header keywords:
+            CX10, CX11, CY10, CY11, IDCSCALE, IDCTHETA, IDCXREF, IDCYREF,
+            IDCV2REF, IDCV3REF
+        """
+        if self.idctab in [None, '', ' ','N/A']:
+            #Keyword idctab is not present in header - check for sip coefficients
+            if header is not None and 'IDCSCALE' in header:
+                self._readModelFromHeader(header)
+            else:
+                print("Distortion model is not available: IDCTAB=None\n")
+                self.idcmodel = None
+        elif not os.path.exists(fileutil.osfn(self.idctab)):
+            if header is not None and 'IDCSCALE' in header:
+                self._readModelFromHeader(header)
+            else:
+                print('Distortion model is not available: IDCTAB file %s not found\n' % self.idctab)
+                self.idcmodel = None
+        else:
+            self.readModelFromIDCTAB(header=header, update=update)
+
+    def _readModelFromHeader(self, header):
+        # Recreate idc model from SIP coefficients and header kw
+        print('Restoring IDC model from SIP coefficients\n')
+        model = models.GeometryModel()
+        cx, cy = coeff_converter.sip2idc(self)
+        model.cx = cx
+        model.cy = cy
+        model.name = "sip"
+        model.norder = header['A_ORDER']
+
+        refpix = {}
+        refpix['XREF'] = header['IDCXREF']
+        refpix['YREF'] = header['IDCYREF']
+        refpix['PSCALE'] = header['IDCSCALE']
+        refpix['V2REF'] = header['IDCV2REF']
+        refpix['V3REF'] = header['IDCV3REF']
+        refpix['THETA'] = header['IDCTHETA']
+        model.refpix = refpix
+
+        self.idcmodel = model
+
+
+    def readModelFromIDCTAB(self, header=None, update=False):
+        """
+        Read distortion model from idc table.
+
+        Parameters
+        ----------
+        header : `astropy.io.fits.Header`
+            fits extension header
+        update : booln (False)
+            if True - save teh following as header keywords:
+            CX10, CX11, CY10, CY11, IDCSCALE, IDCTHETA, IDCXREF, IDCYREF,
+            IDCV2REF, IDCV3REF
+
+        """
+        if self.date_obs == None:
+            print('date_obs not available\n')
+            self.idcmodel = None
+            return
+        if self.filter1 == None and self.filter2 == None:
+            'No filter information available\n'
+            self.idcmodel = None
+            return
+
+        self.idcmodel = models.IDCModel(self.idctab,
+                    chip=self.chip, direction='forward', date=self.date_obs,
+                    filter1=self.filter1, filter2=self.filter2,
+                    offtab=self.offtab, binned=self.binned)
+
+        if self.ltv1 != 0. or self.ltv2 != 0.:
+            self.resetLTV()
+
+        if update:
+            if header==None:
+                print('Update header with IDC model kw requested but header was not provided\n.')
+            else:
+                self._updatehdr(header)
+
+    def resetLTV(self):
+        """
+        Reset LTV values for polarizer data
+
+        The polarizer field is smaller than the detector field.
+        The distortion coefficients are defined for the entire
+        polarizer field and the LTV values are set as with subarray
+        data. This may also be true for other special filters.
+        This is a special case when the observation is considered
+        a subarray in terms of detector field but a full frame in
+        terms of distortion model.
+        To avoid shifting the distortion coefficients the LTV values
+        are reset to 0.
+        """
+        if self.naxis1 == self.idcmodel.refpix['XSIZE'] and  \
+           self.naxis2 == self.idcmodel.refpix['YSIZE']:
+            self.ltv1 = 0.
+            self.ltv2 = 0.
+
+    def wcs2header(self, sip2hdr=False, idc2hdr=True, wcskey=None, relax=False):
+        """
+        Create a `astropy.io.fits.Header` object from WCS keywords.
+
+        If the original header had a CD matrix, return a CD matrix,
+        otherwise return a PC matrix.
+
+        Parameters
+        ----------
+        sip2hdr : bool
+            If True - include SIP coefficients
+        """
+
+        h = self.to_header(key=wcskey, relax=relax)
+        if not wcskey:
+            wcskey = self.wcs.alt
+        if self.wcs.has_cd():
+            h = altwcs.pc2cd(h, key=wcskey)
+
+        if 'wcsname' not in h:
+            if self.idctab is not None:
+                wname = build_default_wcsname(self.idctab)
+            else:
+                wname = 'DEFAULT'
+            h['wcsname{0}'.format(wcskey)] = wname
+
+        if idc2hdr:
+            for card in self._idc2hdr():
+                h[card.keyword + wcskey] = (card.value, card.comment)
+        try:
+            del h['RESTFRQ']
+            del h['RESTWAV']
+        except KeyError: pass
+
+        if sip2hdr and self.sip:
+            for card in self._sip2hdr('a'):
+                h[card.keyword] = (card.value, card.comment)
+            for card in self._sip2hdr('b'):
+                h[card.keyword] = (card.value, card.comment)
+
+            try:
+                ap = self.sip.ap
+            except AssertionError:
+                ap = None
+            try:
+                bp = self.sip.bp
+            except AssertionError:
+                bp = None
+
+            if ap:
+                for card in self._sip2hdr('ap'):
+                    h[card.keyword] = (card.value, card.comment)
+            if bp:
+                for card in self._sip2hdr('bp'):
+                    h[card.keyword] = (card.value, card.comment)
+        return h
+
+    def _sip2hdr(self, k):
+        """
+        Get a set of SIP coefficients in the form of an array
+        and turn them into a `astropy.io.fits.Cardlist`.
+        k - one of 'a', 'b', 'ap', 'bp'
+        """
+
+        cards = [] #fits.CardList()
+        korder = self.sip.__getattribute__(k+'_order')
+        cards.append(fits.Card(keyword=k.upper()+'_ORDER', value=korder))
+        coeffs = self.sip.__getattribute__(k)
+        ind = coeffs.nonzero()
+        for i in range(len(ind[0])):
+            card = fits.Card(keyword=k.upper()+'_'+str(ind[0][i])+'_'+str(ind[1][i]),
+                             value=coeffs[ind[0][i], ind[1][i]])
+            cards.append(card)
+        return cards
+
+    def _idc2hdr(self):
+        # save some of the idc coefficients
+        coeffs = ['ocx10', 'ocx11', 'ocy10', 'ocy11', 'idcscale']
+        cards = [] #fits.CardList()
+        for c in coeffs:
+            try:
+                val = self.__getattribute__(c)
+            except AttributeError:
+                continue
+            if val:
+                cards.append(fits.Card(keyword=c, value=val))
+        return cards
+
+    def pc2cd(self):
+        if not self.wcs.has_pc():
+            self.wcs.pc = self.wcs.get_pc()
+        self.wcs.cd = self.wcs.pc * self.wcs.cdelt[1]
+
+    def all_world2pix(self, *args, **kwargs):
+        """
+        all_world2pix(*arg, accuracy=1.0e-4, maxiter=20, adaptive=False, \
+detect_divergence=True, quiet=False)
+
+        Performs full inverse transformation using iterative solution
+        on full forward transformation with complete distortion model.
+
+        Parameters
+        ----------
+        accuracy : float, optional (Default = 1.0e-4)
+            Required accuracy of the solution. Iteration terminates when the
+            correction to the solution found during the previous iteration
+            is smaller (in the sence of the L2 norm) than `accuracy`\ .
+
+        maxiter : int, optional (Default = 20)
+            Maximum number of iterations allowed to reach the solution.
+
+        adaptive : bool, optional (Default = False)
+            Specifies whether to adaptively select only points that did not
+            converge to a solution whithin the required accuracy for the
+            next iteration. Default is recommended for HST as well as most
+            other instruments.
+
+            .. note::
+               The :py:meth:`all_world2pix` uses a vectorized implementation
+               of the method of consecutive approximations (see `Notes`
+               section below) in which it iterates over *all* input poits
+               *regardless* until the required accuracy has been reached for
+               *all* input points. In some cases it may be possible that
+               *almost all* points have reached the required accuracy but
+               there are only a few of input data points left for which
+               additional iterations may be needed (this depends mostly on the
+               characteristics of the geometric distortions for a given
+               instrument). In this situation it may be
+               advantageous to set `adaptive` = `True`\ in which
+               case :py:meth:`all_world2pix` will continue iterating *only* over
+               the points that have not yet converged to the required
+               accuracy. However, for the HST's ACS/WFC detector, which has
+               the strongest distortions of all HST instruments, testing has
+               shown that enabling this option would lead to a about 10-30\%
+               penalty in computational time (depending on specifics of the
+               image, geometric distortions, and number of input points to be
+               converted). Therefore, for HST instruments,
+               it is recommended to set `adaptive` = `False`\ . The only
+               danger in getting this setting wrong will be a performance
+               penalty.
+
+            .. note::
+               When `detect_divergence` is `True`\ , :py:meth:`all_world2pix` \
+               will automatically switch to the adaptive algorithm once
+               divergence has been detected.
+
+        detect_divergence : bool, optional (Default = True)
+            Specifies whether to perform a more detailed analysis of the
+            convergence to a solution. Normally :py:meth:`all_world2pix`
+            may not achieve the required accuracy
+            if either the `tolerance` or `maxiter` arguments are too low.
+            However, it may happen that for some geometric distortions
+            the conditions of convergence for the the method of consecutive
+            approximations used by :py:meth:`all_world2pix` may not be
+            satisfied, in which case consecutive approximations to the
+            solution will diverge regardless of the `tolerance` or `maxiter`
+            settings.
+
+            When `detect_divergence` is `False`\ , these divergent points
+            will be detected as not having achieved the required accuracy
+            (without further details). In addition, if `adaptive` is `False`
+            then the algorithm will not know that the solution (for specific
+            points) is diverging and will continue iterating and trying to
+            "improve" diverging solutions. This may result in NaN or Inf
+            values in the return results (in addition to a performance
+            penalties). Even when `detect_divergence` is
+            `False`\ , :py:meth:`all_world2pix`\ , at the end of the iterative
+            process, will identify invalid results (NaN or Inf) as "diverging"
+            solutions and will raise :py:class:`NoConvergence` unless
+            the `quiet` parameter is set to `True`\ .
+
+            When `detect_divergence` is `True`\ , :py:meth:`all_world2pix` will
+            detect points for
+            which current correction to the coordinates is larger than
+            the correction applied during the previous iteration **if** the
+            requested accuracy **has not yet been achieved**\ . In this case,
+            if `adaptive` is `True`, these points will be excluded from
+            further iterations and if `adaptive`
+            is `False`\ , :py:meth:`all_world2pix` will automatically
+            switch to the adaptive algorithm.
+
+            .. note::
+               When accuracy has been achieved, small increases in
+               current corrections may be possible due to rounding errors
+               (when `adaptive` is `False`\ ) and such increases
+               will be ignored.
+
+            .. note::
+               Setting `detect_divergence` to `True` will incurr about 5-10\%
+               performance penalty (in our testing on ACS/WFC images).
+               Because the benefits of enabling this feature outweigh
+               the small performance penalty, it is recommended to set
+               `detect_divergence` to `True`\ , unless extensive testing
+               of the distortion models for images from specific
+               instruments show a good stability of the numerical method
+               for a wide range of coordinates (even outside the image
+               itself).
+
+            .. note::
+               Indices of the diverging inverse solutions will be reported
+               in the `divergent` attribute of the
+               raised :py:class:`NoConvergence` object.
+
+        quiet : bool, optional (Default = False)
+            Do not throw :py:class:`NoConvergence` exceptions when the method
+            does not converge to a solution with the required accuracy
+            within a specified number of maximum iterations set by `maxiter`
+            parameter. Instead, simply return the found solution.
+
+        Raises
+        ------
+        NoConvergence
+            The method does not converge to a
+            solution with the required accuracy within a specified number
+            of maximum iterations set by the `maxiter` parameter.
+
+        Notes
+        -----
+        Inputs can either be (RA, Dec, origin) or (RADec, origin) where RA
+        and Dec are 1-D arrays/lists of coordinates and RADec is an
+        array/list of pairs of coordinates.
+
+        Using the method of consecutive approximations we iterate starting
+        with the initial approximation, which is computed using the
+        non-distorion-aware :py:meth:`wcs_world2pix` (or equivalent).
+
+        The :py:meth:`all_world2pix` function uses a vectorized implementation
+        of the method of consecutive approximations and therefore it is
+        highly efficient (>30x) when *all* data points that need to be
+        converted from sky coordinates to image coordinates are passed at
+        *once*\ . Therefore, it is advisable, whenever possible, to pass
+        as input a long array of all points that need to be converted
+        to :py:meth:`all_world2pix` instead of calling :py:meth:`all_world2pix`
+        for each data point. Also see the note to the `adaptive` parameter.
+
+        Examples
+        --------
+        >>> import stwcs
+        >>> from astropy.io import fits
+        >>> hdulist = fits.open('j94f05bgq_flt.fits')
+        >>> w = stwcs.wcsutil.HSTWCS(hdulist, ext=('sci',1))
+        >>> hdulist.close()
+
+        >>> ra, dec = w.all_pix2world([1,2,3],[1,1,1],1); print(ra); print(dec)
+        [ 5.52645241  5.52649277  5.52653313]
+        [-72.05171776 -72.05171295 -72.05170814]
+        >>> radec = w.all_pix2world([[1,1],[2,1],[3,1]],1); print(radec)
+        [[  5.52645241 -72.05171776]
+         [  5.52649277 -72.05171295]
+         [  5.52653313 -72.05170814]]
+        >>> x, y = w.all_world2pix(ra,dec,1)
+        >>> print(x)
+        [ 1.00000233  2.00000232  3.00000233]
+        >>> print(y)
+        [ 0.99999997  0.99999997  0.99999998]
+        >>> xy = w.all_world2pix(radec,1)
+        >>> print(xy)
+        [[ 1.00000233  0.99999997]
+         [ 2.00000232  0.99999997]
+         [ 3.00000233  0.99999998]]
+        >>> xy = w.all_world2pix(radec,1, maxiter=3, accuracy=1.0e-10, \
+quiet=False)
+        NoConvergence: 'HSTWCS.all_world2pix' failed to converge to requested \
+accuracy after 3 iterations.
+
+        >>>
+        Now try to use some diverging data:
+        >>> divradec = w.all_pix2world([[1.0,1.0],[10000.0,50000.0],\
+[3.0,1.0]],1); print(divradec)
+        [[  5.52645241 -72.05171776]
+         [  7.15979392 -70.81405561]
+         [  5.52653313 -72.05170814]]
+
+        >>> try:
+        >>>   xy = w.all_world2pix(divradec,1, maxiter=20, accuracy=1.0e-4, \
+adaptive=False, detect_divergence=True, quiet=False)
+        >>> except stwcs.wcsutil.hstwcs.NoConvergence as e:
+        >>>   print("Indices of diverging points: {}".format(e.divergent))
+        >>>   print("Indices of poorly converging points: {}".format(e.failed2converge))
+        >>>   print("Best solution: {}".format(e.best_solution))
+        >>>   print("Achieved accuracy: {}".format(e.accuracy))
+        >>>   raise e
+        Indices of diverging points:
+        [1]
+        Indices of poorly converging points:
+        None
+        Best solution:
+        [[  1.00006219e+00   9.99999288e-01]
+         [ -1.99440907e+06   1.44308548e+06]
+         [  3.00006257e+00   9.99999316e-01]]
+        Achieved accuracy:
+        [[  5.98554253e-05   6.79918148e-07]
+         [  8.59514088e+11   6.61703754e+11]
+         [  6.02334592e-05   6.59713067e-07]]
+        Traceback (innermost last):
+          File "<console>", line 8, in <module>
+        NoConvergence: 'HSTWCS.all_world2pix' failed to converge to the requested accuracy.
+        After 5 iterations, the solution is diverging at least for one input point.
+
+        >>> try:
+        >>>   xy = w.all_world2pix(divradec,1, maxiter=20, accuracy=1.0e-4, \
+adaptive=False, detect_divergence=False, quiet=False)
+        >>> except stwcs.wcsutil.hstwcs.NoConvergence as e:
+        >>>   print("Indices of diverging points: {}".format(e.divergent))
+        >>>   print("Indices of poorly converging points: {}".format(e.failed2converge))
+        >>>   print("Best solution: {}".format(e.best_solution))
+        >>>   print("Achieved accuracy: {}".format(e.accuracy))
+        >>>   raise e
+        Indices of diverging points:
+        [1]
+        Indices of poorly converging points:
+        None
+        Best solution:
+        [[  1.   1.]
+         [ nan  nan]
+         [  3.   1.]]
+        Achieved accuracy:
+        [[  0.   0.]
+         [ nan  nan]
+         [  0.   0.]]
+        Traceback (innermost last):
+          File "<console>", line 8, in <module>
+        NoConvergence: 'HSTWCS.all_world2pix' failed to converge to the requested accuracy.
+        After 20 iterations, the solution is diverging at least for one input point.
+
+        """
+        #####################################################################
+        ##                     PROCESS ARGUMENTS:                          ##
+        #####################################################################
+        nargs = len(args)
+
+        if nargs == 3:
+            try:
+                ra     = np.asarray(args[0], dtype=np.float64)
+                dec    = np.asarray(args[1], dtype=np.float64)
+                #assert( len(ra.shape) == 1 and len(dec.shape) == 1 )
+                origin = int(args[2])
+                vect1D = True
+            except:
+                raise TypeError("When providing three arguments, they must " \
+                            "be (Ra, Dec, origin) where Ra and Dec are " \
+                            "Nx1 vectors.")
+        elif nargs == 2:
+            try:
+                rd  = np.asarray(args[0], dtype=np.float64)
+                #assert( rd.shape[1] == 2 )
+                ra  = rd[:,0]
+                dec = rd[:,1]
+                origin = int(args[1])
+                vect1D = False
+            except:
+                raise TypeError("When providing two arguments, they must " \
+                            "be (RaDec, origin) where RaDec is a Nx2 array.")
+        else:
+            raise TypeError("Expected 2 or 3 arguments, {:d} given." \
+                            .format(nargs))
+
+        # process optional arguments:
+        accuracy          = kwargs.pop('accuracy', 1.0e-4)
+        maxiter           = kwargs.pop('maxiter', 20)
+        adaptive          = kwargs.pop('adaptive', False)
+        detect_divergence = kwargs.pop('detect_divergence', True)
+        quiet             = kwargs.pop('quiet', False)
+
+        #####################################################################
+        ##                INITIALIZE ITERATIVE PROCESS:                    ##
+        #####################################################################
+        x0, y0 = self.wcs_world2pix(ra, dec, origin) # <-- initial approximation
+                                                  #     (WCS based only)
+
+        # see if an iterative solution is required (when any of the
+        # non-CD-matrix corrections are present). If not required
+        # return initial approximation (x0, y0).
+        if self.sip is None and \
+           self.cpdis1 is None and self.cpdis2 is None and \
+           self.det2im1 is None and self.det2im2 is None:
+            # no non-WCS corrections are detected - return
+            # initial approximation
+            if vect1D:
+                return [x0, y0]
+            else:
+                return np.dstack([x0,y0])[0]
+
+        x  = x0.copy() # 0-order solution
+        y  = y0.copy() # 0-order solution
+
+        # initial correction:
+        dx, dy = self.pix2foc(x, y, origin)
+        # If pix2foc does not apply all the required distortion
+        # corrections then replace the above line with:
+        #r0, d0 = self.all_pix2world(x, y, origin)
+        #dx, dy = self.wcs_world2pix(r0, d0, origin )
+        dx -= x0
+        dy -= y0
+
+        # update initial solution:
+        x -= dx
+        y -= dy
+
+        # norn (L2) squared of the correction:
+        dn2prev   = dx**2+dy**2
+        dn2       = dn2prev
+
+        # prepare for iterative process
+        iterlist  = list(range(1, maxiter+1))
+        accuracy2 = accuracy**2
+        ind       = None
+        inddiv    = None
+
+        npts      = x.shape[0]
+
+        # turn off numpy runtime warnings for 'invalid' and 'over':
+        old_invalid = np.geterr()['invalid']
+        old_over    = np.geterr()['over']
+        np.seterr(invalid = 'ignore', over = 'ignore')
+
+        #####################################################################
+        ##                     NON-ADAPTIVE ITERATIONS:                    ##
+        #####################################################################
+        if not adaptive:
+            for k in iterlist:
+                # check convergence:
+                if np.max(dn2) < accuracy2:
+                    break
+
+                # find correction to the previous solution:
+                dx, dy = self.pix2foc(x, y, origin)
+                # If pix2foc does not apply all the required distortion
+                # corrections then replace the above line with:
+                #r0, d0 = self.all_pix2world(x, y, origin)
+                #dx, dy = self.wcs_world2pix(r0, d0, origin )
+                dx -= x0
+                dy -= y0
+
+                # update norn (L2) squared of the correction:
+                dn2 = dx**2+dy**2
+
+                # check for divergence (we do this in two stages
+                # to optimize performance for the most common
+                # scenario when succesive approximations converge):
+                if detect_divergence:
+                    ind, = np.where(dn2 <= dn2prev)
+                    if ind.shape[0] < npts:
+                        inddiv, = np.where(
+                            np.logical_and(dn2 > dn2prev, dn2 >= accuracy2))
+                        if inddiv.shape[0] > 0:
+                            # apply correction only to the converging points:
+                            x[ind] -= dx[ind]
+                            y[ind] -= dy[ind]
+                            # switch to adaptive iterations:
+                            ind, = np.where((dn2 >= accuracy2) & \
+                                        (dn2 <= dn2prev) & np.isfinite(dn2))
+                            iterlist = iterlist[k:]
+                            adaptive = True
+                            break
+                    #dn2prev[ind] = dn2[ind]
+                    dn2prev = dn2
+
+                # apply correction:
+                x -= dx
+                y -= dy
+
+        #####################################################################
+        ##                      ADAPTIVE ITERATIONS:                       ##
+        #####################################################################
+        if adaptive:
+            if ind is None:
+                ind = np.asarray(list(range(npts)), dtype=np.int64)
+
+            for k in iterlist:
+                # check convergence:
+                if ind.shape[0] == 0:
+                    break
+
+                # find correction to the previous solution:
+                dx[ind], dy[ind] = self.pix2foc(x[ind], y[ind], origin)
+                # If pix2foc does not apply all the required distortion
+                # corrections then replace the above line with:
+                #r0[ind], d0[ind] = self.all_pix2world(x[ind], y[ind], origin)
+                #dx[ind], dy[ind] = self.wcs_world2pix(r0[ind], d0[ind], origin)
+                dx[ind] -= x0[ind]
+                dy[ind] -= y0[ind]
+
+                # update norn (L2) squared of the correction:
+                dn2 = dx**2+dy**2
+
+                # update indices of elements that still need correction:
+                if detect_divergence:
+                    ind, = np.where((dn2 >= accuracy2) & (dn2 <= dn2prev))
+                    #ind = ind[np.where((dn2[ind] >= accuracy2) & (dn2[ind] <= dn2prev))]
+                    dn2prev[ind] = dn2[ind]
+                else:
+                    ind, = np.where(dn2 >= accuracy2)
+                    #ind = ind[np.where(dn2[ind] >= accuracy2)]
+
+                # apply correction:
+                x[ind] -= dx[ind]
+                y[ind] -= dy[ind]
+
+        #####################################################################
+        ##         FINAL DETECTION OF INVALID, DIVERGING,                  ##
+        ##         AND FAILED-TO-CONVERGE POINTS                           ##
+        #####################################################################
+        # Identify diverging and/or invalid points:
+        invalid = (((~np.isfinite(y)) | (~np.isfinite(x)) | \
+                    (~np.isfinite(dn2))) & \
+                   (np.isfinite(ra)) &  (np.isfinite(dec)))
+        # When detect_divergence==False, dn2prev is outdated (it is the
+        # norm^2 of the very first correction). Still better than nothing...
+        inddiv, = np.where(((dn2 >= accuracy2) & (dn2 > dn2prev)) | invalid)
+        if inddiv.shape[0] == 0:
+            inddiv = None
+        # identify points that did not converge within
+        # 'maxiter' iterations:
+        if k >= maxiter:
+            ind,= np.where((dn2 >= accuracy2) & (dn2 <= dn2prev) & (~invalid))
+            if ind.shape[0] == 0:
+                ind = None
+        else:
+            ind = None
+
+        #####################################################################
+        ##      RAISE EXCEPTION IF DIVERGING OR TOO SLOWLY CONVERGING      ##
+        ##      DATA POINTS HAVE BEEN DETECTED:                            ##
+        #####################################################################
+        # raise exception if diverging or too slowly converging
+        if (ind is not None or inddiv is not None) and not quiet:
+            if vect1D:
+                sol  = [x, y]
+                err  = [np.abs(dx), np.abs(dy)]
+            else:
+                sol  = np.dstack( [x, y] )[0]
+                err  = np.dstack( [np.abs(dx), np.abs(dy)] )[0]
+
+            # restore previous numpy error settings:
+            np.seterr(invalid = old_invalid, over = old_over)
+
+            if inddiv is None:
+                raise NoConvergence("'HSTWCS.all_world2pix' failed to "        \
+                    "converge to the requested accuracy after {:d} "         \
+                    "iterations.".format(k), best_solution = sol,            \
+                    accuracy = err, niter = k, failed2converge = ind,        \
+                    divergent = None)
+            else:
+                raise NoConvergence("'HSTWCS.all_world2pix' failed to "        \
+                    "converge to the requested accuracy.{0:s}"               \
+                    "After {1:d} iterations, the solution is diverging "     \
+                    "at least for one input point."                          \
+                    .format(os.linesep, k), best_solution = sol,             \
+                    accuracy = err, niter = k, failed2converge = ind,        \
+                    divergent = inddiv)
+
+        #####################################################################
+        ##             FINALIZE AND FORMAT DATA FOR RETURN:                ##
+        #####################################################################
+        # restore previous numpy error settings:
+        np.seterr(invalid = old_invalid, over = old_over)
+
+        if vect1D:
+            return [x, y]
+        else:
+            return np.dstack( [x, y] )[0]
+
+    def _updatehdr(self, ext_hdr):
+        #kw2add : OCX10, OCX11, OCY10, OCY11
+        # record the model in the header for use by pydrizzle
+        ext_hdr['OCX10'] = self.idcmodel.cx[1,0]
+        ext_hdr['OCX11'] = self.idcmodel.cx[1,1]
+        ext_hdr['OCY10'] = self.idcmodel.cy[1,0]
+        ext_hdr['OCY11'] = self.idcmodel.cy[1,1]
+        ext_hdr['IDCSCALE'] = self.idcmodel.refpix['PSCALE']
+        ext_hdr['IDCTHETA'] = self.idcmodel.refpix['THETA']
+        ext_hdr['IDCXREF'] = self.idcmodel.refpix['XREF']
+        ext_hdr['IDCYREF'] = self.idcmodel.refpix['YREF']
+        ext_hdr['IDCV2REF'] =  self.idcmodel.refpix['V2REF']
+        ext_hdr['IDCV3REF'] = self.idcmodel.refpix['V3REF']
+
+    def printwcs(self):
+        """
+        Print the basic WCS keywords.
+        """
+        print('WCS Keywords\n')
+        print('CD_11  CD_12: %r %r' % (self.wcs.cd[0,0],  self.wcs.cd[0,1]))
+        print('CD_21  CD_22: %r %r' % (self.wcs.cd[1,0],  self.wcs.cd[1,1]))
+        print('CRVAL    : %r %r' % (self.wcs.crval[0], self.wcs.crval[1]))
+        print('CRPIX    : %r %r' % (self.wcs.crpix[0], self.wcs.crpix[1]))
+        print('NAXIS    : %d %d' % (self.naxis1, self.naxis2))
+        print('Plate Scale : %r' % self.pscale)
+        print('ORIENTAT : %r' % self.orientat)
+
+
+def determine_refframe(phdr):
+    """
+    Determine the reference frame in standard FITS WCS terms.
+
+    Parameters
+    ----------
+    phdr : `astropy.io.fits.Header`
+        Primary Header of an HST observation
+
+    In HST images the reference frame is recorded in the primary extension as REFFRAME.
+    Values are "GSC1" which means FK5 or ICRS (for GSC2 observations).
+    """
+    try:
+        refframe = phdr['REFFRAME']
+    except KeyError:
+        refframe = " "
+    if refframe == "GSC1":
+        refframe = "FK5"
+    return refframe