from __future__ import division # confidence high import os #from pywcs import WCS from astropy.wcs import WCS from astropy.io import fits as pyfits #import pyfits import instruments from stwcs.distortion import models, coeff_converter import altwcs import numpy as np from stsci.tools import fileutil import getinput import mappings 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_sky2pix`\ . 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`. slow_conv : 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 points (i.e., if the required accuracy has been achieved for all points) then `slow_conv` 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.slow_conv = kwargs.pop('slow_conv', 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: string or PyFITS HDUList object or None a file name, e.g j9irw4b1q_flt.fits a fully qualified filename[EXTNAME,EXTNUM], e.g. j9irw4b1q_flt.fits[sci,1] a pyfits file object, e.g pyfits.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 in ['IRAF/ARTDATA','',' ','N/A']: self.instrument = 'DEFAULT' else: self.instrument = instrument_name WCS.__init__(self, ehdr, fobj=phdu, minerr=self.minerr, key=self.wcskey) # If input was a pyfits HDUList object, it's the user's # responsibility to close it, otherwise, it's closed here. if not isinstance(fobj, pyfits.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() def readIDCCoeffs(self, header): """ Reads in first order IDCTAB coefficients if present in the header """ coeffs = ['ocx10', 'ocx11', 'ocy10', 'ocy11', 'idcscale'] 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: pyfits.Header primary header ext_hdr: pyfits.Header extension header """ if self.instrument in 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: pyfits.Header fits extension header update: boolean (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: pyfits.Header fits extension header update: boolean (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 pyfits.Header object from WCS keywords. If the original header had a CD matrix, return a CD matrix, otherwise return a PC matrix. Parameters ---------- sip2hdr: boolean 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.update('wcsname'+wcskey, value=wname) if idc2hdr: for card in self._idc2hdr(): h.update(card.key+wcskey, value=card.value, comment=card.comment) try: del h['RESTFRQ'] del h['RESTWAV'] except KeyError: pass if sip2hdr and self.sip: for card in self._sip2hdr('a'): h.update(card.key,value=card.value,comment=card.comment) for card in self._sip2hdr('b'): h.update(card.key,value=card.value,comment=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.update(card.key,value=card.value,comment=card.comment) if bp: for card in self._sip2hdr('bp'): h.update(card.key,value=card.value,comment=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 pyfits.Cardlist. k - one of 'a', 'b', 'ap', 'bp' """ cards = pyfits.CardList() korder = self.sip.__getattribute__(k+'_order') cards.append(pyfits.Card(key=k.upper()+'_ORDER', value=korder)) coeffs = self.sip.__getattribute__(k) ind = coeffs.nonzero() for i in range(len(ind[0])): card = pyfits.Card(key=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 = pyfits.CardList() for c in coeffs: try: val = self.__getattribute__(c) except AttributeError: continue if val: cards.append(pyfits.Card(key=c, value=val)) return cards def pc2cd(self): self.wcs.cd = self.wcs.pc.copy() def all_sky2pix(self, *args, **kwargs): """ all_sky2pix(*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_sky2pix` 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_sky2pix` 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_sky2pix` \ 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_sky2pix` 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_sky2pix` 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_sky2pix`\ , 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_sky2pix` 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_sky2pix` 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` exception 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_sky2pix` (or equivalent). The :py:meth:`all_sky2pix` 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_sky2pix` instead of calling :py:meth:`all_sky2pix` for each data point. Also see the note to the `adaptive` parameter. See Also -------- A detailed description of the algorithm is available on astropy's GitHub page, issue `#2373 `_\ . Examples -------- >>> import stwcs, pyfits >>> hdulist = pyfits.open('j94f05bgq_flt.fits') >>> w = stwcs.wcsutil.HSTWCS(hdulist, ext=('sci',1)) >>> hdulist.close() >>> ra, dec = w.all_pix2sky([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_pix2sky([[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_sky2pix(ra,dec,1) >>> print(x) [ 1.00000233 2.00000232 3.00000233] >>> print(y) [ 0.99999997 0.99999997 0.99999998] >>> xy = w.all_sky2pix(radec,1) >>> print(xy) [[ 1.00000233 0.99999997] [ 2.00000232 0.99999997] [ 3.00000233 0.99999998]] >>> xy = w.all_sky2pix(radec,1, maxiter=3, accuracy=1.0e-10, \ quiet=False) NoConvergence: 'HSTWCS.all_sky2pix' failed to converge to requested \ accuracy after 3 iterations. >>> Now try to use some diverging data: >>> divradec = w.all_pix2sky([[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_sky2pix(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.slow_conv)) >>> 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 "", line 8, in NoConvergence: 'HSTWCS.all_sky2pix' failed to converge to the \ requested accuracy. After 5 iterations, the solution is diverging at least for one \ input point. >>> try: >>> xy = w.all_sky2pix(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.slow_conv)) >>> 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 "", line 8, in NoConvergence: 'HSTWCS.all_sky2pix' 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_sky2pix(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_pix2sky(x, y, origin) #dx, dy = self.wcs_sky2pix(r0, d0, origin ) dx -= x0 dy -= y0 # update initial solution: x -= dx y -= dy # norm (L2) squared of the correction: dn2prev = dx**2+dy**2 dn2 = dn2prev # prepare for iterative process iterlist = 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_pix2sky(x, y, origin) #dx, dy = self.wcs_sky2pix(r0, d0, origin ) dx -= x0 dy -= y0 # update norm (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: if np.any(dn2 > dn2prev): 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(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_pix2sky(x[ind], y[ind], origin) #dx[ind], dy[ind] = self.wcs_sky2pix(r0[ind], d0[ind], origin) dx[ind] -= x0[ind] dy[ind] -= y0[ind] # update norm (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_sky2pix' failed to " \ "converge to the requested accuracy after {:d} " \ "iterations.".format(k), best_solution = sol, \ accuracy = err, niter = k, slow_conv = ind, \ divergent = None) else: raise NoConvergence("'HSTWCS.all_sky2pix' 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, slow_conv = 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.update('OCX10', self.idcmodel.cx[1,0]) ext_hdr.update('OCX11', self.idcmodel.cx[1,1]) ext_hdr.update('OCY10', self.idcmodel.cy[1,0]) ext_hdr.update('OCY11', self.idcmodel.cy[1,1]) ext_hdr.update('IDCSCALE', self.idcmodel.refpix['PSCALE']) ext_hdr.update('IDCTHETA', self.idcmodel.refpix['THETA']) ext_hdr.update('IDCXREF', self.idcmodel.refpix['XREF']) ext_hdr.update('IDCYREF', self.idcmodel.refpix['YREF']) ext_hdr.update('IDCV2REF', self.idcmodel.refpix['V2REF']) ext_hdr.update('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