1# Licensed under a 3-clause BSD style license - see LICENSE.rst 2 3# Under the hood, there are 3 separate classes that perform different 4# parts of the transformation: 5# 6# - `~astropy.wcs.Wcsprm`: Is a direct wrapper of the core WCS 7# functionality in `wcslib`_. (This includes TPV and TPD 8# polynomial distortion, but not SIP distortion). 9# 10# - `~astropy.wcs.Sip`: Handles polynomial distortion as defined in the 11# `SIP`_ convention. 12# 13# - `~astropy.wcs.DistortionLookupTable`: Handles `distortion paper`_ 14# lookup tables. 15# 16# Additionally, the class `WCS` aggregates all of these transformations 17# together in a pipeline: 18# 19# - Detector to image plane correction (by a pair of 20# `~astropy.wcs.DistortionLookupTable` objects). 21# 22# - `SIP`_ distortion correction (by an underlying `~astropy.wcs.Sip` 23# object) 24# 25# - `distortion paper`_ table-lookup correction (by a pair of 26# `~astropy.wcs.DistortionLookupTable` objects). 27# 28# - `wcslib`_ WCS transformation (by a `~astropy.wcs.Wcsprm` object) 29 30# STDLIB 31import copy 32import uuid 33import io 34import itertools 35import os 36import re 37import textwrap 38import warnings 39import builtins 40 41# THIRD-PARTY 42import numpy as np 43 44# LOCAL 45from astropy import log 46from astropy.io import fits 47from . import docstrings 48from . import _wcs 49 50from astropy import units as u 51from astropy.utils.compat import possible_filename 52from astropy.utils.exceptions import AstropyWarning, AstropyUserWarning, AstropyDeprecationWarning 53from astropy.utils.decorators import deprecated_renamed_argument 54 55# Mix-in class that provides the APE 14 API 56from .wcsapi.fitswcs import FITSWCSAPIMixin, SlicedFITSWCS 57 58__all__ = ['FITSFixedWarning', 'WCS', 'find_all_wcs', 59 'DistortionLookupTable', 'Sip', 'Tabprm', 'Wcsprm', 'Auxprm', 60 'Wtbarr', 'WCSBase', 'validate', 'WcsError', 'SingularMatrixError', 61 'InconsistentAxisTypesError', 'InvalidTransformError', 62 'InvalidCoordinateError', 'NoSolutionError', 63 'InvalidSubimageSpecificationError', 'NoConvergence', 64 'NonseparableSubimageCoordinateSystemError', 65 'NoWcsKeywordsFoundError', 'InvalidTabularParametersError'] 66 67 68__doctest_skip__ = ['WCS.all_world2pix'] 69 70 71if _wcs is not None: 72 _parsed_version = _wcs.__version__.split('.') 73 if int(_parsed_version[0]) == 5 and int(_parsed_version[1]) < 8: 74 raise ImportError( 75 "astropy.wcs is built with wcslib {0}, but only versions 5.8 and " 76 "later on the 5.x series are known to work. The version of wcslib " 77 "that ships with astropy may be used.") 78 79 if not _wcs._sanity_check(): 80 raise RuntimeError( 81 "astropy.wcs did not pass its sanity check for your build " 82 "on your platform.") 83 84 WCSBase = _wcs._Wcs 85 DistortionLookupTable = _wcs.DistortionLookupTable 86 Sip = _wcs.Sip 87 Wcsprm = _wcs.Wcsprm 88 Auxprm = _wcs.Auxprm 89 Tabprm = _wcs.Tabprm 90 Wtbarr = _wcs.Wtbarr 91 WcsError = _wcs.WcsError 92 SingularMatrixError = _wcs.SingularMatrixError 93 InconsistentAxisTypesError = _wcs.InconsistentAxisTypesError 94 InvalidTransformError = _wcs.InvalidTransformError 95 InvalidCoordinateError = _wcs.InvalidCoordinateError 96 NoSolutionError = _wcs.NoSolutionError 97 InvalidSubimageSpecificationError = _wcs.InvalidSubimageSpecificationError 98 NonseparableSubimageCoordinateSystemError = _wcs.NonseparableSubimageCoordinateSystemError 99 NoWcsKeywordsFoundError = _wcs.NoWcsKeywordsFoundError 100 InvalidTabularParametersError = _wcs.InvalidTabularParametersError 101 102 # Copy all the constants from the C extension into this module's namespace 103 for key, val in _wcs.__dict__.items(): 104 if key.startswith(('WCSSUB_', 'WCSHDR_', 'WCSHDO_', 'WCSCOMPARE_')): 105 locals()[key] = val 106 __all__.append(key) 107 108 # Set coordinate extraction callback for WCS -TAB: 109 def _load_tab_bintable(hdulist, extnam, extver, extlev, kind, ttype, row, ndim): 110 arr = hdulist[(extnam, extver)].data[ttype][row - 1] 111 112 if arr.ndim != ndim: 113 if kind == 'c' and ndim == 2: 114 arr = arr.reshape((arr.size, 1)) 115 else: 116 raise ValueError("Bad TDIM") 117 118 return np.ascontiguousarray(arr, dtype=np.double) 119 120 _wcs.set_wtbarr_fitsio_callback(_load_tab_bintable) 121 122else: 123 WCSBase = object 124 Wcsprm = object 125 DistortionLookupTable = object 126 Sip = object 127 Tabprm = object 128 Wtbarr = object 129 WcsError = None 130 SingularMatrixError = None 131 InconsistentAxisTypesError = None 132 InvalidTransformError = None 133 InvalidCoordinateError = None 134 NoSolutionError = None 135 InvalidSubimageSpecificationError = None 136 NonseparableSubimageCoordinateSystemError = None 137 NoWcsKeywordsFoundError = None 138 InvalidTabularParametersError = None 139 140 141# Additional relax bit flags 142WCSHDO_SIP = 0x80000 143 144# Regular expression defining SIP keyword It matches keyword that starts with A 145# or B, optionally followed by P, followed by an underscore then a number in 146# range of 0-19, followed by an underscore and another number in range of 0-19. 147# Keyword optionally ends with a capital letter. 148SIP_KW = re.compile('''^[AB]P?_1?[0-9]_1?[0-9][A-Z]?$''') 149 150 151def _parse_keysel(keysel): 152 keysel_flags = 0 153 if keysel is not None: 154 for element in keysel: 155 if element.lower() == 'image': 156 keysel_flags |= _wcs.WCSHDR_IMGHEAD 157 elif element.lower() == 'binary': 158 keysel_flags |= _wcs.WCSHDR_BIMGARR 159 elif element.lower() == 'pixel': 160 keysel_flags |= _wcs.WCSHDR_PIXLIST 161 else: 162 raise ValueError( 163 "keysel must be a list of 'image', 'binary' " + 164 "and/or 'pixel'") 165 else: 166 keysel_flags = -1 167 168 return keysel_flags 169 170 171class NoConvergence(Exception): 172 """ 173 An error class used to report non-convergence and/or divergence 174 of numerical methods. It is used to report errors in the 175 iterative solution used by 176 the :py:meth:`~astropy.wcs.WCS.all_world2pix`. 177 178 Attributes 179 ---------- 180 181 best_solution : `numpy.ndarray` 182 Best solution achieved by the numerical method. 183 184 accuracy : `numpy.ndarray` 185 Accuracy of the ``best_solution``. 186 187 niter : `int` 188 Number of iterations performed by the numerical method 189 to compute ``best_solution``. 190 191 divergent : None, `numpy.ndarray` 192 Indices of the points in ``best_solution`` array 193 for which the solution appears to be divergent. If the 194 solution does not diverge, ``divergent`` will be set to `None`. 195 196 slow_conv : None, `numpy.ndarray` 197 Indices of the solutions in ``best_solution`` array 198 for which the solution failed to converge within the 199 specified maximum number of iterations. If there are no 200 non-converging solutions (i.e., if the required accuracy 201 has been achieved for all input data points) 202 then ``slow_conv`` will be set to `None`. 203 204 """ 205 206 def __init__(self, *args, best_solution=None, accuracy=None, niter=None, 207 divergent=None, slow_conv=None, **kwargs): 208 super().__init__(*args) 209 210 self.best_solution = best_solution 211 self.accuracy = accuracy 212 self.niter = niter 213 self.divergent = divergent 214 self.slow_conv = slow_conv 215 216 if kwargs: 217 warnings.warn("Function received unexpected arguments ({}) these " 218 "are ignored but will raise an Exception in the " 219 "future.".format(list(kwargs)), 220 AstropyDeprecationWarning) 221 222 223class FITSFixedWarning(AstropyWarning): 224 """ 225 The warning raised when the contents of the FITS header have been 226 modified to be standards compliant. 227 """ 228 pass 229 230 231class WCS(FITSWCSAPIMixin, WCSBase): 232 """WCS objects perform standard WCS transformations, and correct for 233 `SIP`_ and `distortion paper`_ table-lookup transformations, based 234 on the WCS keywords and supplementary data read from a FITS file. 235 236 See also: https://docs.astropy.org/en/stable/wcs/ 237 238 Parameters 239 ---------- 240 header : `~astropy.io.fits.Header`, `~astropy.io.fits.hdu.image.PrimaryHDU`, `~astropy.io.fits.hdu.image.ImageHDU`, str, dict-like, or None, optional 241 If *header* is not provided or None, the object will be 242 initialized to default values. 243 244 fobj : `~astropy.io.fits.HDUList`, optional 245 It is needed when header keywords point to a `distortion 246 paper`_ lookup table stored in a different extension. 247 248 key : str, optional 249 The name of a particular WCS transform to use. This may be 250 either ``' '`` or ``'A'``-``'Z'`` and corresponds to the 251 ``\"a\"`` part of the ``CTYPEia`` cards. *key* may only be 252 provided if *header* is also provided. 253 254 minerr : float, optional 255 The minimum value a distortion correction must have in order 256 to be applied. If the value of ``CQERRja`` is smaller than 257 *minerr*, the corresponding distortion is not applied. 258 259 relax : bool or int, optional 260 Degree of permissiveness: 261 262 - `True` (default): Admit all recognized informal extensions 263 of the WCS standard. 264 265 - `False`: Recognize only FITS keywords defined by the 266 published WCS standard. 267 268 - `int`: a bit field selecting specific extensions to accept. 269 See :ref:`astropy:relaxread` for details. 270 271 naxis : int or sequence, optional 272 Extracts specific coordinate axes using 273 :meth:`~astropy.wcs.Wcsprm.sub`. If a header is provided, and 274 *naxis* is not ``None``, *naxis* will be passed to 275 :meth:`~astropy.wcs.Wcsprm.sub` in order to select specific 276 axes from the header. See :meth:`~astropy.wcs.Wcsprm.sub` for 277 more details about this parameter. 278 279 keysel : sequence of str, optional 280 A sequence of flags used to select the keyword types 281 considered by wcslib. When ``None``, only the standard image 282 header keywords are considered (and the underlying wcspih() C 283 function is called). To use binary table image array or pixel 284 list keywords, *keysel* must be set. 285 286 Each element in the list should be one of the following 287 strings: 288 289 - 'image': Image header keywords 290 291 - 'binary': Binary table image array keywords 292 293 - 'pixel': Pixel list keywords 294 295 Keywords such as ``EQUIna`` or ``RFRQna`` that are common to 296 binary table image arrays and pixel lists (including 297 ``WCSNna`` and ``TWCSna``) are selected by both 'binary' and 298 'pixel'. 299 300 colsel : sequence of int, optional 301 A sequence of table column numbers used to restrict the WCS 302 transformations considered to only those pertaining to the 303 specified columns. If `None`, there is no restriction. 304 305 fix : bool, optional 306 When `True` (default), call `~astropy.wcs.Wcsprm.fix` on 307 the resulting object to fix any non-standard uses in the 308 header. `FITSFixedWarning` Warnings will be emitted if any 309 changes were made. 310 311 translate_units : str, optional 312 Specify which potentially unsafe translations of non-standard 313 unit strings to perform. By default, performs none. See 314 `WCS.fix` for more information about this parameter. Only 315 effective when ``fix`` is `True`. 316 317 Raises 318 ------ 319 MemoryError 320 Memory allocation failed. 321 322 ValueError 323 Invalid key. 324 325 KeyError 326 Key not found in FITS header. 327 328 ValueError 329 Lookup table distortion present in the header but *fobj* was 330 not provided. 331 332 Notes 333 ----- 334 335 1. astropy.wcs supports arbitrary *n* dimensions for the core WCS 336 (the transformations handled by WCSLIB). However, the 337 `distortion paper`_ lookup table and `SIP`_ distortions must be 338 two dimensional. Therefore, if you try to create a WCS object 339 where the core WCS has a different number of dimensions than 2 340 and that object also contains a `distortion paper`_ lookup 341 table or `SIP`_ distortion, a `ValueError` 342 exception will be raised. To avoid this, consider using the 343 *naxis* kwarg to select two dimensions from the core WCS. 344 345 2. The number of coordinate axes in the transformation is not 346 determined directly from the ``NAXIS`` keyword but instead from 347 the highest of: 348 349 - ``NAXIS`` keyword 350 351 - ``WCSAXESa`` keyword 352 353 - The highest axis number in any parameterized WCS keyword. 354 The keyvalue, as well as the keyword, must be 355 syntactically valid otherwise it will not be considered. 356 357 If none of these keyword types is present, i.e. if the header 358 only contains auxiliary WCS keywords for a particular 359 coordinate representation, then no coordinate description is 360 constructed for it. 361 362 The number of axes, which is set as the ``naxis`` member, may 363 differ for different coordinate representations of the same 364 image. 365 366 3. When the header includes duplicate keywords, in most cases the 367 last encountered is used. 368 369 4. `~astropy.wcs.Wcsprm.set` is called immediately after 370 construction, so any invalid keywords or transformations will 371 be raised by the constructor, not when subsequently calling a 372 transformation method. 373 374 """ # noqa: E501 375 376 def __init__(self, header=None, fobj=None, key=' ', minerr=0.0, 377 relax=True, naxis=None, keysel=None, colsel=None, 378 fix=True, translate_units='', _do_set=True): 379 close_fds = [] 380 381 if header is None: 382 if naxis is None: 383 naxis = 2 384 wcsprm = _wcs.Wcsprm(header=None, key=key, 385 relax=relax, naxis=naxis) 386 self.naxis = wcsprm.naxis 387 # Set some reasonable defaults. 388 det2im = (None, None) 389 cpdis = (None, None) 390 sip = None 391 else: 392 keysel_flags = _parse_keysel(keysel) 393 394 if isinstance(header, (str, bytes)): 395 try: 396 is_path = (possible_filename(header) and 397 os.path.exists(header)) 398 except (OSError, ValueError): 399 is_path = False 400 401 if is_path: 402 if fobj is not None: 403 raise ValueError( 404 "Can not provide both a FITS filename to " 405 "argument 1 and a FITS file object to argument 2") 406 fobj = fits.open(header) 407 close_fds.append(fobj) 408 header = fobj[0].header 409 elif isinstance(header, fits.hdu.image._ImageBaseHDU): 410 header = header.header 411 elif not isinstance(header, fits.Header): 412 try: 413 # Accept any dict-like object 414 orig_header = header 415 header = fits.Header() 416 for dict_key in orig_header.keys(): 417 header[dict_key] = orig_header[dict_key] 418 except TypeError: 419 raise TypeError( 420 "header must be a string, an astropy.io.fits.Header " 421 "object, or a dict-like object") 422 423 if isinstance(header, fits.Header): 424 header_string = header.tostring().rstrip() 425 else: 426 header_string = header 427 428 # Importantly, header is a *copy* of the passed-in header 429 # because we will be modifying it 430 if isinstance(header_string, str): 431 header_bytes = header_string.encode('ascii') 432 header_string = header_string 433 else: 434 header_bytes = header_string 435 header_string = header_string.decode('ascii') 436 437 if not (fobj is None or isinstance(fobj, fits.HDUList)): 438 raise AssertionError("'fobj' must be either None or an " 439 "astropy.io.fits.HDUList object.") 440 441 est_naxis = 2 442 try: 443 tmp_header = fits.Header.fromstring(header_string) 444 self._remove_sip_kw(tmp_header) 445 tmp_header_bytes = tmp_header.tostring().rstrip() 446 if isinstance(tmp_header_bytes, str): 447 tmp_header_bytes = tmp_header_bytes.encode('ascii') 448 tmp_wcsprm = _wcs.Wcsprm(header=tmp_header_bytes, key=key, 449 relax=relax, keysel=keysel_flags, 450 colsel=colsel, warnings=False, 451 hdulist=fobj) 452 if naxis is not None: 453 try: 454 tmp_wcsprm = tmp_wcsprm.sub(naxis) 455 except ValueError: 456 pass 457 est_naxis = tmp_wcsprm.naxis if tmp_wcsprm.naxis else 2 458 459 except _wcs.NoWcsKeywordsFoundError: 460 pass 461 462 self.naxis = est_naxis 463 464 header = fits.Header.fromstring(header_string) 465 466 det2im = self._read_det2im_kw(header, fobj, err=minerr) 467 cpdis = self._read_distortion_kw( 468 header, fobj, dist='CPDIS', err=minerr) 469 sip = self._read_sip_kw(header, wcskey=key) 470 self._remove_sip_kw(header) 471 472 header_string = header.tostring() 473 header_string = header_string.replace('END' + ' ' * 77, '') 474 475 if isinstance(header_string, str): 476 header_bytes = header_string.encode('ascii') 477 header_string = header_string 478 else: 479 header_bytes = header_string 480 header_string = header_string.decode('ascii') 481 482 try: 483 wcsprm = _wcs.Wcsprm(header=header_bytes, key=key, 484 relax=relax, keysel=keysel_flags, 485 colsel=colsel, hdulist=fobj) 486 except _wcs.NoWcsKeywordsFoundError: 487 # The header may have SIP or distortions, but no core 488 # WCS. That isn't an error -- we want a "default" 489 # (identity) core Wcs transformation in that case. 490 if colsel is None: 491 wcsprm = _wcs.Wcsprm(header=None, key=key, 492 relax=relax, keysel=keysel_flags, 493 colsel=colsel, hdulist=fobj) 494 else: 495 raise 496 497 if naxis is not None: 498 wcsprm = wcsprm.sub(naxis) 499 self.naxis = wcsprm.naxis 500 501 if (wcsprm.naxis != 2 and 502 (det2im[0] or det2im[1] or cpdis[0] or cpdis[1] or sip)): 503 raise ValueError( 504 """ 505FITS WCS distortion paper lookup tables and SIP distortions only work 506in 2 dimensions. However, WCSLIB has detected {} dimensions in the 507core WCS keywords. To use core WCS in conjunction with FITS WCS 508distortion paper lookup tables or SIP distortion, you must select or 509reduce these to 2 dimensions using the naxis kwarg. 510""".format(wcsprm.naxis)) 511 512 header_naxis = header.get('NAXIS', None) 513 if header_naxis is not None and header_naxis < wcsprm.naxis: 514 warnings.warn( 515 "The WCS transformation has more axes ({:d}) than the " 516 "image it is associated with ({:d})".format( 517 wcsprm.naxis, header_naxis), FITSFixedWarning) 518 519 self._get_naxis(header) 520 WCSBase.__init__(self, sip, cpdis, wcsprm, det2im) 521 522 if fix: 523 if header is None: 524 with warnings.catch_warnings(): 525 warnings.simplefilter('ignore', FITSFixedWarning) 526 self.fix(translate_units=translate_units) 527 else: 528 self.fix(translate_units=translate_units) 529 530 if _do_set: 531 self.wcs.set() 532 533 for fd in close_fds: 534 fd.close() 535 536 self._pixel_bounds = None 537 538 def __copy__(self): 539 new_copy = self.__class__() 540 WCSBase.__init__(new_copy, self.sip, 541 (self.cpdis1, self.cpdis2), 542 self.wcs, 543 (self.det2im1, self.det2im2)) 544 new_copy.__dict__.update(self.__dict__) 545 return new_copy 546 547 def __deepcopy__(self, memo): 548 from copy import deepcopy 549 550 new_copy = self.__class__() 551 new_copy.naxis = deepcopy(self.naxis, memo) 552 WCSBase.__init__(new_copy, deepcopy(self.sip, memo), 553 (deepcopy(self.cpdis1, memo), 554 deepcopy(self.cpdis2, memo)), 555 deepcopy(self.wcs, memo), 556 (deepcopy(self.det2im1, memo), 557 deepcopy(self.det2im2, memo))) 558 for key, val in self.__dict__.items(): 559 new_copy.__dict__[key] = deepcopy(val, memo) 560 return new_copy 561 562 def copy(self): 563 """ 564 Return a shallow copy of the object. 565 566 Convenience method so user doesn't have to import the 567 :mod:`copy` stdlib module. 568 569 .. warning:: 570 Use `deepcopy` instead of `copy` unless you know why you need a 571 shallow copy. 572 """ 573 return copy.copy(self) 574 575 def deepcopy(self): 576 """ 577 Return a deep copy of the object. 578 579 Convenience method so user doesn't have to import the 580 :mod:`copy` stdlib module. 581 """ 582 return copy.deepcopy(self) 583 584 def sub(self, axes=None): 585 586 copy = self.deepcopy() 587 588 # We need to know which axes have been dropped, but there is no easy 589 # way to do this with the .sub function, so instead we assign UUIDs to 590 # the CNAME parameters in copy.wcs. We can later access the original 591 # CNAME properties from self.wcs. 592 cname_uuid = [str(uuid.uuid4()) for i in range(copy.wcs.naxis)] 593 copy.wcs.cname = cname_uuid 594 595 # Subset the WCS 596 copy.wcs = copy.wcs.sub(axes) 597 copy.naxis = copy.wcs.naxis 598 599 # Construct a list of dimensions from the original WCS in the order 600 # in which they appear in the final WCS. 601 keep = [cname_uuid.index(cname) if cname in cname_uuid else None 602 for cname in copy.wcs.cname] 603 604 # Restore the original CNAMEs 605 copy.wcs.cname = ['' if i is None else self.wcs.cname[i] for i in keep] 606 607 # Subset pixel_shape and pixel_bounds 608 if self.pixel_shape: 609 copy.pixel_shape = tuple([None if i is None else self.pixel_shape[i] for i in keep]) 610 if self.pixel_bounds: 611 copy.pixel_bounds = [None if i is None else self.pixel_bounds[i] for i in keep] 612 613 return copy 614 615 if _wcs is not None: 616 sub.__doc__ = _wcs.Wcsprm.sub.__doc__ 617 618 def _fix_scamp(self): 619 """ 620 Remove SCAMP's PVi_m distortion parameters if SIP distortion parameters 621 are also present. Some projects (e.g., Palomar Transient Factory) 622 convert SCAMP's distortion parameters (which abuse the PVi_m cards) to 623 SIP. However, wcslib gets confused by the presence of both SCAMP and 624 SIP distortion parameters. 625 626 See https://github.com/astropy/astropy/issues/299. 627 """ 628 # Nothing to be done if no WCS attached 629 if self.wcs is None: 630 return 631 632 # Nothing to be done if no PV parameters attached 633 pv = self.wcs.get_pv() 634 if not pv: 635 return 636 637 # Nothing to be done if axes don't use SIP distortion parameters 638 if self.sip is None: 639 return 640 641 # Nothing to be done if any radial terms are present... 642 # Loop over list to find any radial terms. 643 # Certain values of the `j' index are used for storing 644 # radial terms; refer to Equation (1) in 645 # <http://web.ipac.caltech.edu/staff/shupe/reprints/SIP_to_PV_SPIE2012.pdf>. 646 pv = np.asarray(pv) 647 # Loop over distinct values of `i' index 648 for i in set(pv[:, 0]): 649 # Get all values of `j' index for this value of `i' index 650 js = set(pv[:, 1][pv[:, 0] == i]) 651 # Find max value of `j' index 652 max_j = max(js) 653 for j in (3, 11, 23, 39): 654 if j < max_j and j in js: 655 return 656 657 self.wcs.set_pv([]) 658 warnings.warn("Removed redundant SCAMP distortion parameters " + 659 "because SIP parameters are also present", FITSFixedWarning) 660 661 def fix(self, translate_units='', naxis=None): 662 """ 663 Perform the fix operations from wcslib, and warn about any 664 changes it has made. 665 666 Parameters 667 ---------- 668 translate_units : str, optional 669 Specify which potentially unsafe translations of 670 non-standard unit strings to perform. By default, 671 performs none. 672 673 Although ``"S"`` is commonly used to represent seconds, 674 its translation to ``"s"`` is potentially unsafe since the 675 standard recognizes ``"S"`` formally as Siemens, however 676 rarely that may be used. The same applies to ``"H"`` for 677 hours (Henry), and ``"D"`` for days (Debye). 678 679 This string controls what to do in such cases, and is 680 case-insensitive. 681 682 - If the string contains ``"s"``, translate ``"S"`` to 683 ``"s"``. 684 685 - If the string contains ``"h"``, translate ``"H"`` to 686 ``"h"``. 687 688 - If the string contains ``"d"``, translate ``"D"`` to 689 ``"d"``. 690 691 Thus ``''`` doesn't do any unsafe translations, whereas 692 ``'shd'`` does all of them. 693 694 naxis : int array, optional 695 Image axis lengths. If this array is set to zero or 696 ``None``, then `~astropy.wcs.Wcsprm.cylfix` will not be 697 invoked. 698 """ 699 if self.wcs is not None: 700 self._fix_scamp() 701 fixes = self.wcs.fix(translate_units, naxis) 702 for key, val in fixes.items(): 703 if val != "No change": 704 if (key == 'datfix' and '1858-11-17' in val and 705 not np.count_nonzero(self.wcs.mjdref)): 706 continue 707 warnings.warn( 708 ("'{0}' made the change '{1}'."). 709 format(key, val), 710 FITSFixedWarning) 711 712 def calc_footprint(self, header=None, undistort=True, axes=None, center=True): 713 """ 714 Calculates the footprint of the image on the sky. 715 716 A footprint is defined as the positions of the corners of the 717 image on the sky after all available distortions have been 718 applied. 719 720 Parameters 721 ---------- 722 header : `~astropy.io.fits.Header` object, optional 723 Used to get ``NAXIS1`` and ``NAXIS2`` 724 header and axes are mutually exclusive, alternative ways 725 to provide the same information. 726 727 undistort : bool, optional 728 If `True`, take SIP and distortion lookup table into 729 account 730 731 axes : (int, int), optional 732 If provided, use the given sequence as the shape of the 733 image. Otherwise, use the ``NAXIS1`` and ``NAXIS2`` 734 keywords from the header that was used to create this 735 `WCS` object. 736 737 center : bool, optional 738 If `True` use the center of the pixel, otherwise use the corner. 739 740 Returns 741 ------- 742 coord : (4, 2) array of (*x*, *y*) coordinates. 743 The order is clockwise starting with the bottom left corner. 744 """ 745 if axes is not None: 746 naxis1, naxis2 = axes 747 else: 748 if header is None: 749 try: 750 # classes that inherit from WCS and define naxis1/2 751 # do not require a header parameter 752 naxis1, naxis2 = self.pixel_shape 753 except (AttributeError, TypeError): 754 warnings.warn( 755 "Need a valid header in order to calculate footprint\n", AstropyUserWarning) 756 return None 757 else: 758 naxis1 = header.get('NAXIS1', None) 759 naxis2 = header.get('NAXIS2', None) 760 761 if naxis1 is None or naxis2 is None: 762 raise ValueError( 763 "Image size could not be determined.") 764 765 if center: 766 corners = np.array([[1, 1], 767 [1, naxis2], 768 [naxis1, naxis2], 769 [naxis1, 1]], dtype=np.float64) 770 else: 771 corners = np.array([[0.5, 0.5], 772 [0.5, naxis2 + 0.5], 773 [naxis1 + 0.5, naxis2 + 0.5], 774 [naxis1 + 0.5, 0.5]], dtype=np.float64) 775 776 if undistort: 777 return self.all_pix2world(corners, 1) 778 else: 779 return self.wcs_pix2world(corners, 1) 780 781 def _read_det2im_kw(self, header, fobj, err=0.0): 782 """ 783 Create a `distortion paper`_ type lookup table for detector to 784 image plane correction. 785 """ 786 if fobj is None: 787 return (None, None) 788 789 if not isinstance(fobj, fits.HDUList): 790 return (None, None) 791 792 try: 793 axiscorr = header['AXISCORR'] 794 d2imdis = self._read_d2im_old_format(header, fobj, axiscorr) 795 return d2imdis 796 except KeyError: 797 pass 798 799 dist = 'D2IMDIS' 800 d_kw = 'D2IM' 801 err_kw = 'D2IMERR' 802 tables = {} 803 for i in range(1, self.naxis + 1): 804 d_error = header.get(err_kw + str(i), 0.0) 805 if d_error < err: 806 tables[i] = None 807 continue 808 distortion = dist + str(i) 809 if distortion in header: 810 dis = header[distortion].lower() 811 if dis == 'lookup': 812 del header[distortion] 813 assert isinstance(fobj, fits.HDUList), ( 814 'An astropy.io.fits.HDUList' 815 'is required for Lookup table distortion.') 816 dp = (d_kw + str(i)).strip() 817 dp_extver_key = dp + '.EXTVER' 818 if dp_extver_key in header: 819 d_extver = header[dp_extver_key] 820 del header[dp_extver_key] 821 else: 822 d_extver = 1 823 dp_axis_key = dp + f'.AXIS.{i:d}' 824 if i == header[dp_axis_key]: 825 d_data = fobj['D2IMARR', d_extver].data 826 else: 827 d_data = (fobj['D2IMARR', d_extver].data).transpose() 828 del header[dp_axis_key] 829 d_header = fobj['D2IMARR', d_extver].header 830 d_crpix = (d_header.get('CRPIX1', 0.0), d_header.get('CRPIX2', 0.0)) 831 d_crval = (d_header.get('CRVAL1', 0.0), d_header.get('CRVAL2', 0.0)) 832 d_cdelt = (d_header.get('CDELT1', 1.0), d_header.get('CDELT2', 1.0)) 833 d_lookup = DistortionLookupTable(d_data, d_crpix, 834 d_crval, d_cdelt) 835 tables[i] = d_lookup 836 else: 837 warnings.warn('Polynomial distortion is not implemented.\n', AstropyUserWarning) 838 for key in set(header): 839 if key.startswith(dp + '.'): 840 del header[key] 841 else: 842 tables[i] = None 843 if not tables: 844 return (None, None) 845 else: 846 return (tables.get(1), tables.get(2)) 847 848 def _read_d2im_old_format(self, header, fobj, axiscorr): 849 warnings.warn( 850 "The use of ``AXISCORR`` for D2IM correction has been deprecated." 851 "`~astropy.wcs` will read in files with ``AXISCORR`` but ``to_fits()`` will write " 852 "out files without it.", 853 AstropyDeprecationWarning) 854 cpdis = [None, None] 855 crpix = [0., 0.] 856 crval = [0., 0.] 857 cdelt = [1., 1.] 858 try: 859 d2im_data = fobj[('D2IMARR', 1)].data 860 except KeyError: 861 return (None, None) 862 except AttributeError: 863 return (None, None) 864 865 d2im_data = np.array([d2im_data]) 866 d2im_hdr = fobj[('D2IMARR', 1)].header 867 naxis = d2im_hdr['NAXIS'] 868 869 for i in range(1, naxis + 1): 870 crpix[i - 1] = d2im_hdr.get('CRPIX' + str(i), 0.0) 871 crval[i - 1] = d2im_hdr.get('CRVAL' + str(i), 0.0) 872 cdelt[i - 1] = d2im_hdr.get('CDELT' + str(i), 1.0) 873 874 cpdis = DistortionLookupTable(d2im_data, crpix, crval, cdelt) 875 876 if axiscorr == 1: 877 return (cpdis, None) 878 elif axiscorr == 2: 879 return (None, cpdis) 880 else: 881 warnings.warn("Expected AXISCORR to be 1 or 2", AstropyUserWarning) 882 return (None, None) 883 884 def _write_det2im(self, hdulist): 885 """ 886 Writes a `distortion paper`_ type lookup table to the given 887 `~astropy.io.fits.HDUList`. 888 """ 889 890 if self.det2im1 is None and self.det2im2 is None: 891 return 892 dist = 'D2IMDIS' 893 d_kw = 'D2IM' 894 895 def write_d2i(num, det2im): 896 if det2im is None: 897 return 898 899 hdulist[0].header[f'{dist}{num:d}'] = ( 900 'LOOKUP', 'Detector to image correction type') 901 hdulist[0].header[f'{d_kw}{num:d}.EXTVER'] = ( 902 num, 'Version number of WCSDVARR extension') 903 hdulist[0].header[f'{d_kw}{num:d}.NAXES'] = ( 904 len(det2im.data.shape), 'Number of independent variables in D2IM function') 905 906 for i in range(det2im.data.ndim): 907 jth = {1: '1st', 2: '2nd', 3: '3rd'}.get(i + 1, f'{i + 1}th') 908 hdulist[0].header[f'{d_kw}{num:d}.AXIS.{i + 1:d}'] = ( 909 i + 1, f'Axis number of the {jth} variable in a D2IM function') 910 911 image = fits.ImageHDU(det2im.data, name='D2IMARR') 912 header = image.header 913 914 header['CRPIX1'] = (det2im.crpix[0], 915 'Coordinate system reference pixel') 916 header['CRPIX2'] = (det2im.crpix[1], 917 'Coordinate system reference pixel') 918 header['CRVAL1'] = (det2im.crval[0], 919 'Coordinate system value at reference pixel') 920 header['CRVAL2'] = (det2im.crval[1], 921 'Coordinate system value at reference pixel') 922 header['CDELT1'] = (det2im.cdelt[0], 923 'Coordinate increment along axis') 924 header['CDELT2'] = (det2im.cdelt[1], 925 'Coordinate increment along axis') 926 image.ver = int(hdulist[0].header[f'{d_kw}{num:d}.EXTVER']) 927 hdulist.append(image) 928 write_d2i(1, self.det2im1) 929 write_d2i(2, self.det2im2) 930 931 def _read_distortion_kw(self, header, fobj, dist='CPDIS', err=0.0): 932 """ 933 Reads `distortion paper`_ table-lookup keywords and data, and 934 returns a 2-tuple of `~astropy.wcs.DistortionLookupTable` 935 objects. 936 937 If no `distortion paper`_ keywords are found, ``(None, None)`` 938 is returned. 939 """ 940 if isinstance(header, (str, bytes)): 941 return (None, None) 942 943 if dist == 'CPDIS': 944 d_kw = 'DP' 945 err_kw = 'CPERR' 946 else: 947 d_kw = 'DQ' 948 err_kw = 'CQERR' 949 950 tables = {} 951 for i in range(1, self.naxis + 1): 952 d_error_key = err_kw + str(i) 953 if d_error_key in header: 954 d_error = header[d_error_key] 955 del header[d_error_key] 956 else: 957 d_error = 0.0 958 if d_error < err: 959 tables[i] = None 960 continue 961 distortion = dist + str(i) 962 if distortion in header: 963 dis = header[distortion].lower() 964 del header[distortion] 965 if dis == 'lookup': 966 if not isinstance(fobj, fits.HDUList): 967 raise ValueError('an astropy.io.fits.HDUList is ' 968 'required for Lookup table distortion.') 969 dp = (d_kw + str(i)).strip() 970 dp_extver_key = dp + '.EXTVER' 971 if dp_extver_key in header: 972 d_extver = header[dp_extver_key] 973 del header[dp_extver_key] 974 else: 975 d_extver = 1 976 dp_axis_key = dp + f'.AXIS.{i:d}' 977 if i == header[dp_axis_key]: 978 d_data = fobj['WCSDVARR', d_extver].data 979 else: 980 d_data = (fobj['WCSDVARR', d_extver].data).transpose() 981 del header[dp_axis_key] 982 d_header = fobj['WCSDVARR', d_extver].header 983 d_crpix = (d_header.get('CRPIX1', 0.0), 984 d_header.get('CRPIX2', 0.0)) 985 d_crval = (d_header.get('CRVAL1', 0.0), 986 d_header.get('CRVAL2', 0.0)) 987 d_cdelt = (d_header.get('CDELT1', 1.0), 988 d_header.get('CDELT2', 1.0)) 989 d_lookup = DistortionLookupTable(d_data, d_crpix, d_crval, d_cdelt) 990 tables[i] = d_lookup 991 992 for key in set(header): 993 if key.startswith(dp + '.'): 994 del header[key] 995 else: 996 warnings.warn('Polynomial distortion is not implemented.\n', AstropyUserWarning) 997 else: 998 tables[i] = None 999 1000 if not tables: 1001 return (None, None) 1002 else: 1003 return (tables.get(1), tables.get(2)) 1004 1005 def _write_distortion_kw(self, hdulist, dist='CPDIS'): 1006 """ 1007 Write out `distortion paper`_ keywords to the given 1008 `~astropy.io.fits.HDUList`. 1009 """ 1010 if self.cpdis1 is None and self.cpdis2 is None: 1011 return 1012 1013 if dist == 'CPDIS': 1014 d_kw = 'DP' 1015 else: 1016 d_kw = 'DQ' 1017 1018 def write_dist(num, cpdis): 1019 if cpdis is None: 1020 return 1021 1022 hdulist[0].header[f'{dist}{num:d}'] = ( 1023 'LOOKUP', 'Prior distortion function type') 1024 hdulist[0].header[f'{d_kw}{num:d}.EXTVER'] = ( 1025 num, 'Version number of WCSDVARR extension') 1026 hdulist[0].header[f'{d_kw}{num:d}.NAXES'] = ( 1027 len(cpdis.data.shape), f'Number of independent variables in {dist} function') 1028 1029 for i in range(cpdis.data.ndim): 1030 jth = {1: '1st', 2: '2nd', 3: '3rd'}.get(i + 1, f'{i + 1}th') 1031 hdulist[0].header[f'{d_kw}{num:d}.AXIS.{i + 1:d}'] = ( 1032 i + 1, 1033 f'Axis number of the {jth} variable in a {dist} function') 1034 1035 image = fits.ImageHDU(cpdis.data, name='WCSDVARR') 1036 header = image.header 1037 1038 header['CRPIX1'] = (cpdis.crpix[0], 'Coordinate system reference pixel') 1039 header['CRPIX2'] = (cpdis.crpix[1], 'Coordinate system reference pixel') 1040 header['CRVAL1'] = (cpdis.crval[0], 'Coordinate system value at reference pixel') 1041 header['CRVAL2'] = (cpdis.crval[1], 'Coordinate system value at reference pixel') 1042 header['CDELT1'] = (cpdis.cdelt[0], 'Coordinate increment along axis') 1043 header['CDELT2'] = (cpdis.cdelt[1], 'Coordinate increment along axis') 1044 image.ver = int(hdulist[0].header[f'{d_kw}{num:d}.EXTVER']) 1045 hdulist.append(image) 1046 1047 write_dist(1, self.cpdis1) 1048 write_dist(2, self.cpdis2) 1049 1050 def _remove_sip_kw(self, header): 1051 """ 1052 Remove SIP information from a header. 1053 """ 1054 # Never pass SIP coefficients to wcslib 1055 # CTYPE must be passed with -SIP to wcslib 1056 for key in set(m.group() for m in map(SIP_KW.match, list(header)) 1057 if m is not None): 1058 del header[key] 1059 1060 def _read_sip_kw(self, header, wcskey=""): 1061 """ 1062 Reads `SIP`_ header keywords and returns a `~astropy.wcs.Sip` 1063 object. 1064 1065 If no `SIP`_ header keywords are found, ``None`` is returned. 1066 """ 1067 if isinstance(header, (str, bytes)): 1068 # TODO: Parse SIP from a string without pyfits around 1069 return None 1070 1071 if "A_ORDER" in header and header['A_ORDER'] > 1: 1072 if "B_ORDER" not in header: 1073 raise ValueError( 1074 "A_ORDER provided without corresponding B_ORDER " 1075 "keyword for SIP distortion") 1076 1077 m = int(header["A_ORDER"]) 1078 a = np.zeros((m + 1, m + 1), np.double) 1079 for i in range(m + 1): 1080 for j in range(m - i + 1): 1081 key = f"A_{i}_{j}" 1082 if key in header: 1083 a[i, j] = header[key] 1084 del header[key] 1085 1086 m = int(header["B_ORDER"]) 1087 if m > 1: 1088 b = np.zeros((m + 1, m + 1), np.double) 1089 for i in range(m + 1): 1090 for j in range(m - i + 1): 1091 key = f"B_{i}_{j}" 1092 if key in header: 1093 b[i, j] = header[key] 1094 del header[key] 1095 else: 1096 a = None 1097 b = None 1098 1099 del header['A_ORDER'] 1100 del header['B_ORDER'] 1101 1102 ctype = [header[f'CTYPE{nax}{wcskey}'] for nax in range(1, self.naxis + 1)] 1103 if any(not ctyp.endswith('-SIP') for ctyp in ctype): 1104 message = """ 1105 Inconsistent SIP distortion information is present in the FITS header and the WCS object: 1106 SIP coefficients were detected, but CTYPE is missing a "-SIP" suffix. 1107 astropy.wcs is using the SIP distortion coefficients, 1108 therefore the coordinates calculated here might be incorrect. 1109 1110 If you do not want to apply the SIP distortion coefficients, 1111 please remove the SIP coefficients from the FITS header or the 1112 WCS object. As an example, if the image is already distortion-corrected 1113 (e.g., drizzled) then distortion components should not apply and the SIP 1114 coefficients should be removed. 1115 1116 While the SIP distortion coefficients are being applied here, if that was indeed the intent, 1117 for consistency please append "-SIP" to the CTYPE in the FITS header or the WCS object. 1118 1119 """ # noqa: E501 1120 log.info(message) 1121 elif "B_ORDER" in header and header['B_ORDER'] > 1: 1122 raise ValueError( 1123 "B_ORDER provided without corresponding A_ORDER " + 1124 "keyword for SIP distortion") 1125 else: 1126 a = None 1127 b = None 1128 1129 if "AP_ORDER" in header and header['AP_ORDER'] > 1: 1130 if "BP_ORDER" not in header: 1131 raise ValueError( 1132 "AP_ORDER provided without corresponding BP_ORDER " 1133 "keyword for SIP distortion") 1134 1135 m = int(header["AP_ORDER"]) 1136 ap = np.zeros((m + 1, m + 1), np.double) 1137 for i in range(m + 1): 1138 for j in range(m - i + 1): 1139 key = f"AP_{i}_{j}" 1140 if key in header: 1141 ap[i, j] = header[key] 1142 del header[key] 1143 1144 m = int(header["BP_ORDER"]) 1145 if m > 1: 1146 bp = np.zeros((m + 1, m + 1), np.double) 1147 for i in range(m + 1): 1148 for j in range(m - i + 1): 1149 key = f"BP_{i}_{j}" 1150 if key in header: 1151 bp[i, j] = header[key] 1152 del header[key] 1153 else: 1154 ap = None 1155 bp = None 1156 1157 del header['AP_ORDER'] 1158 del header['BP_ORDER'] 1159 elif "BP_ORDER" in header and header['BP_ORDER'] > 1: 1160 raise ValueError( 1161 "BP_ORDER provided without corresponding AP_ORDER " 1162 "keyword for SIP distortion") 1163 else: 1164 ap = None 1165 bp = None 1166 1167 if a is None and b is None and ap is None and bp is None: 1168 return None 1169 1170 if f"CRPIX1{wcskey}" not in header or f"CRPIX2{wcskey}" not in header: 1171 raise ValueError( 1172 "Header has SIP keywords without CRPIX keywords") 1173 1174 crpix1 = header.get(f"CRPIX1{wcskey}") 1175 crpix2 = header.get(f"CRPIX2{wcskey}") 1176 1177 return Sip(a, b, ap, bp, (crpix1, crpix2)) 1178 1179 def _write_sip_kw(self): 1180 """ 1181 Write out SIP keywords. Returns a dictionary of key-value 1182 pairs. 1183 """ 1184 if self.sip is None: 1185 return {} 1186 1187 keywords = {} 1188 1189 def write_array(name, a): 1190 if a is None: 1191 return 1192 size = a.shape[0] 1193 trdir = 'sky to detector' if name[-1] == 'P' else 'detector to sky' 1194 comment = ('SIP polynomial order, axis {:d}, {:s}' 1195 .format(ord(name[0]) - ord('A'), trdir)) 1196 keywords[f'{name}_ORDER'] = size - 1, comment 1197 1198 comment = 'SIP distortion coefficient' 1199 for i in range(size): 1200 for j in range(size - i): 1201 if a[i, j] != 0.0: 1202 keywords[ 1203 f'{name}_{i:d}_{j:d}'] = a[i, j], comment 1204 1205 write_array('A', self.sip.a) 1206 write_array('B', self.sip.b) 1207 write_array('AP', self.sip.ap) 1208 write_array('BP', self.sip.bp) 1209 1210 return keywords 1211 1212 def _denormalize_sky(self, sky): 1213 if self.wcs.lngtyp != 'RA': 1214 raise ValueError( 1215 "WCS does not have longitude type of 'RA', therefore " + 1216 "(ra, dec) data can not be used as input") 1217 if self.wcs.lattyp != 'DEC': 1218 raise ValueError( 1219 "WCS does not have longitude type of 'DEC', therefore " + 1220 "(ra, dec) data can not be used as input") 1221 if self.wcs.naxis == 2: 1222 if self.wcs.lng == 0 and self.wcs.lat == 1: 1223 return sky 1224 elif self.wcs.lng == 1 and self.wcs.lat == 0: 1225 # Reverse the order of the columns 1226 return sky[:, ::-1] 1227 else: 1228 raise ValueError( 1229 "WCS does not have longitude and latitude celestial " + 1230 "axes, therefore (ra, dec) data can not be used as input") 1231 else: 1232 if self.wcs.lng < 0 or self.wcs.lat < 0: 1233 raise ValueError( 1234 "WCS does not have both longitude and latitude " 1235 "celestial axes, therefore (ra, dec) data can not be " + 1236 "used as input") 1237 out = np.zeros((sky.shape[0], self.wcs.naxis)) 1238 out[:, self.wcs.lng] = sky[:, 0] 1239 out[:, self.wcs.lat] = sky[:, 1] 1240 return out 1241 1242 def _normalize_sky(self, sky): 1243 if self.wcs.lngtyp != 'RA': 1244 raise ValueError( 1245 "WCS does not have longitude type of 'RA', therefore " + 1246 "(ra, dec) data can not be returned") 1247 if self.wcs.lattyp != 'DEC': 1248 raise ValueError( 1249 "WCS does not have longitude type of 'DEC', therefore " + 1250 "(ra, dec) data can not be returned") 1251 if self.wcs.naxis == 2: 1252 if self.wcs.lng == 0 and self.wcs.lat == 1: 1253 return sky 1254 elif self.wcs.lng == 1 and self.wcs.lat == 0: 1255 # Reverse the order of the columns 1256 return sky[:, ::-1] 1257 else: 1258 raise ValueError( 1259 "WCS does not have longitude and latitude celestial " 1260 "axes, therefore (ra, dec) data can not be returned") 1261 else: 1262 if self.wcs.lng < 0 or self.wcs.lat < 0: 1263 raise ValueError( 1264 "WCS does not have both longitude and latitude celestial " 1265 "axes, therefore (ra, dec) data can not be returned") 1266 out = np.empty((sky.shape[0], 2)) 1267 out[:, 0] = sky[:, self.wcs.lng] 1268 out[:, 1] = sky[:, self.wcs.lat] 1269 return out 1270 1271 def _array_converter(self, func, sky, *args, ra_dec_order=False): 1272 """ 1273 A helper function to support reading either a pair of arrays 1274 or a single Nx2 array. 1275 """ 1276 1277 def _return_list_of_arrays(axes, origin): 1278 if any([x.size == 0 for x in axes]): 1279 return axes 1280 1281 try: 1282 axes = np.broadcast_arrays(*axes) 1283 except ValueError: 1284 raise ValueError( 1285 "Coordinate arrays are not broadcastable to each other") 1286 1287 xy = np.hstack([x.reshape((x.size, 1)) for x in axes]) 1288 1289 if ra_dec_order and sky == 'input': 1290 xy = self._denormalize_sky(xy) 1291 output = func(xy, origin) 1292 if ra_dec_order and sky == 'output': 1293 output = self._normalize_sky(output) 1294 return (output[:, 0].reshape(axes[0].shape), 1295 output[:, 1].reshape(axes[0].shape)) 1296 return [output[:, i].reshape(axes[0].shape) 1297 for i in range(output.shape[1])] 1298 1299 def _return_single_array(xy, origin): 1300 if xy.shape[-1] != self.naxis: 1301 raise ValueError( 1302 "When providing two arguments, the array must be " 1303 "of shape (N, {})".format(self.naxis)) 1304 if 0 in xy.shape: 1305 return xy 1306 if ra_dec_order and sky == 'input': 1307 xy = self._denormalize_sky(xy) 1308 result = func(xy, origin) 1309 if ra_dec_order and sky == 'output': 1310 result = self._normalize_sky(result) 1311 return result 1312 1313 if len(args) == 2: 1314 try: 1315 xy, origin = args 1316 xy = np.asarray(xy) 1317 origin = int(origin) 1318 except Exception: 1319 raise TypeError( 1320 "When providing two arguments, they must be " 1321 "(coords[N][{}], origin)".format(self.naxis)) 1322 if xy.shape == () or len(xy.shape) == 1: 1323 return _return_list_of_arrays([xy], origin) 1324 return _return_single_array(xy, origin) 1325 1326 elif len(args) == self.naxis + 1: 1327 axes = args[:-1] 1328 origin = args[-1] 1329 try: 1330 axes = [np.asarray(x) for x in axes] 1331 origin = int(origin) 1332 except Exception: 1333 raise TypeError( 1334 "When providing more than two arguments, they must be " + 1335 "a 1-D array for each axis, followed by an origin.") 1336 1337 return _return_list_of_arrays(axes, origin) 1338 1339 raise TypeError( 1340 "WCS projection has {0} dimensions, so expected 2 (an Nx{0} array " 1341 "and the origin argument) or {1} arguments (the position in each " 1342 "dimension, and the origin argument). Instead, {2} arguments were " 1343 "given.".format( 1344 self.naxis, self.naxis + 1, len(args))) 1345 1346 def all_pix2world(self, *args, **kwargs): 1347 return self._array_converter( 1348 self._all_pix2world, 'output', *args, **kwargs) 1349 all_pix2world.__doc__ = """ 1350 Transforms pixel coordinates to world coordinates. 1351 1352 Performs all of the following in series: 1353 1354 - Detector to image plane correction (if present in the 1355 FITS file) 1356 1357 - `SIP`_ distortion correction (if present in the FITS 1358 file) 1359 1360 - `distortion paper`_ table-lookup correction (if present 1361 in the FITS file) 1362 1363 - `wcslib`_ "core" WCS transformation 1364 1365 Parameters 1366 ---------- 1367 {} 1368 1369 For a transformation that is not two-dimensional, the 1370 two-argument form must be used. 1371 1372 {} 1373 1374 Returns 1375 ------- 1376 1377 {} 1378 1379 Notes 1380 ----- 1381 The order of the axes for the result is determined by the 1382 ``CTYPEia`` keywords in the FITS header, therefore it may not 1383 always be of the form (*ra*, *dec*). The 1384 `~astropy.wcs.Wcsprm.lat`, `~astropy.wcs.Wcsprm.lng`, 1385 `~astropy.wcs.Wcsprm.lattyp` and `~astropy.wcs.Wcsprm.lngtyp` 1386 members can be used to determine the order of the axes. 1387 1388 Raises 1389 ------ 1390 MemoryError 1391 Memory allocation failed. 1392 1393 SingularMatrixError 1394 Linear transformation matrix is singular. 1395 1396 InconsistentAxisTypesError 1397 Inconsistent or unrecognized coordinate axis types. 1398 1399 ValueError 1400 Invalid parameter value. 1401 1402 ValueError 1403 Invalid coordinate transformation parameters. 1404 1405 ValueError 1406 x- and y-coordinate arrays are not the same size. 1407 1408 InvalidTransformError 1409 Invalid coordinate transformation parameters. 1410 1411 InvalidTransformError 1412 Ill-conditioned coordinate transformation parameters. 1413 """.format(docstrings.TWO_OR_MORE_ARGS('naxis', 8), 1414 docstrings.RA_DEC_ORDER(8), 1415 docstrings.RETURNS('sky coordinates, in degrees', 8)) 1416 1417 def wcs_pix2world(self, *args, **kwargs): 1418 if self.wcs is None: 1419 raise ValueError("No basic WCS settings were created.") 1420 return self._array_converter( 1421 lambda xy, o: self.wcs.p2s(xy, o)['world'], 1422 'output', *args, **kwargs) 1423 wcs_pix2world.__doc__ = """ 1424 Transforms pixel coordinates to world coordinates by doing 1425 only the basic `wcslib`_ transformation. 1426 1427 No `SIP`_ or `distortion paper`_ table lookup correction is 1428 applied. To perform distortion correction, see 1429 `~astropy.wcs.WCS.all_pix2world`, 1430 `~astropy.wcs.WCS.sip_pix2foc`, `~astropy.wcs.WCS.p4_pix2foc`, 1431 or `~astropy.wcs.WCS.pix2foc`. 1432 1433 Parameters 1434 ---------- 1435 {} 1436 1437 For a transformation that is not two-dimensional, the 1438 two-argument form must be used. 1439 1440 {} 1441 1442 Returns 1443 ------- 1444 1445 {} 1446 1447 Raises 1448 ------ 1449 MemoryError 1450 Memory allocation failed. 1451 1452 SingularMatrixError 1453 Linear transformation matrix is singular. 1454 1455 InconsistentAxisTypesError 1456 Inconsistent or unrecognized coordinate axis types. 1457 1458 ValueError 1459 Invalid parameter value. 1460 1461 ValueError 1462 Invalid coordinate transformation parameters. 1463 1464 ValueError 1465 x- and y-coordinate arrays are not the same size. 1466 1467 InvalidTransformError 1468 Invalid coordinate transformation parameters. 1469 1470 InvalidTransformError 1471 Ill-conditioned coordinate transformation parameters. 1472 1473 Notes 1474 ----- 1475 The order of the axes for the result is determined by the 1476 ``CTYPEia`` keywords in the FITS header, therefore it may not 1477 always be of the form (*ra*, *dec*). The 1478 `~astropy.wcs.Wcsprm.lat`, `~astropy.wcs.Wcsprm.lng`, 1479 `~astropy.wcs.Wcsprm.lattyp` and `~astropy.wcs.Wcsprm.lngtyp` 1480 members can be used to determine the order of the axes. 1481 1482 """.format(docstrings.TWO_OR_MORE_ARGS('naxis', 8), 1483 docstrings.RA_DEC_ORDER(8), 1484 docstrings.RETURNS('world coordinates, in degrees', 8)) 1485 1486 def _all_world2pix(self, world, origin, tolerance, maxiter, adaptive, 1487 detect_divergence, quiet): 1488 # ############################################################ 1489 # # DESCRIPTION OF THE NUMERICAL METHOD ## 1490 # ############################################################ 1491 # In this section I will outline the method of solving 1492 # the inverse problem of converting world coordinates to 1493 # pixel coordinates (*inverse* of the direct transformation 1494 # `all_pix2world`) and I will summarize some of the aspects 1495 # of the method proposed here and some of the issues of the 1496 # original `all_world2pix` (in relation to this method) 1497 # discussed in https://github.com/astropy/astropy/issues/1977 1498 # A more detailed discussion can be found here: 1499 # https://github.com/astropy/astropy/pull/2373 1500 # 1501 # 1502 # ### Background ### 1503 # 1504 # 1505 # I will refer here to the [SIP Paper] 1506 # (http://fits.gsfc.nasa.gov/registry/sip/SIP_distortion_v1_0.pdf). 1507 # According to this paper, the effect of distortions as 1508 # described in *their* equation (1) is: 1509 # 1510 # (1) x = CD*(u+f(u)), 1511 # 1512 # where `x` is a *vector* of "intermediate spherical 1513 # coordinates" (equivalent to (x,y) in the paper) and `u` 1514 # is a *vector* of "pixel coordinates", and `f` is a vector 1515 # function describing geometrical distortions 1516 # (see equations 2 and 3 in SIP Paper. 1517 # However, I prefer to use `w` for "intermediate world 1518 # coordinates", `x` for pixel coordinates, and assume that 1519 # transformation `W` performs the **linear** 1520 # (CD matrix + projection onto celestial sphere) part of the 1521 # conversion from pixel coordinates to world coordinates. 1522 # Then we can re-write (1) as: 1523 # 1524 # (2) w = W*(x+f(x)) = T(x) 1525 # 1526 # In `astropy.wcs.WCS` transformation `W` is represented by 1527 # the `wcs_pix2world` member, while the combined ("total") 1528 # transformation (linear part + distortions) is performed by 1529 # `all_pix2world`. Below I summarize the notations and their 1530 # equivalents in `astropy.wcs.WCS`: 1531 # 1532 # | Equation term | astropy.WCS/meaning | 1533 # | ------------- | ---------------------------- | 1534 # | `x` | pixel coordinates | 1535 # | `w` | world coordinates | 1536 # | `W` | `wcs_pix2world()` | 1537 # | `W^{-1}` | `wcs_world2pix()` | 1538 # | `T` | `all_pix2world()` | 1539 # | `x+f(x)` | `pix2foc()` | 1540 # 1541 # 1542 # ### Direct Solving of Equation (2) ### 1543 # 1544 # 1545 # In order to find the pixel coordinates that correspond to 1546 # given world coordinates `w`, it is necessary to invert 1547 # equation (2): `x=T^{-1}(w)`, or solve equation `w==T(x)` 1548 # for `x`. However, this approach has the following 1549 # disadvantages: 1550 # 1. It requires unnecessary transformations (see next 1551 # section). 1552 # 2. It is prone to "RA wrapping" issues as described in 1553 # https://github.com/astropy/astropy/issues/1977 1554 # (essentially because `all_pix2world` may return points with 1555 # a different phase than user's input `w`). 1556 # 1557 # 1558 # ### Description of the Method Used here ### 1559 # 1560 # 1561 # By applying inverse linear WCS transformation (`W^{-1}`) 1562 # to both sides of equation (2) and introducing notation `x'` 1563 # (prime) for the pixels coordinates obtained from the world 1564 # coordinates by applying inverse *linear* WCS transformation 1565 # ("focal plane coordinates"): 1566 # 1567 # (3) x' = W^{-1}(w) 1568 # 1569 # we obtain the following equation: 1570 # 1571 # (4) x' = x+f(x), 1572 # 1573 # or, 1574 # 1575 # (5) x = x'-f(x) 1576 # 1577 # This equation is well suited for solving using the method 1578 # of fixed-point iterations 1579 # (http://en.wikipedia.org/wiki/Fixed-point_iteration): 1580 # 1581 # (6) x_{i+1} = x'-f(x_i) 1582 # 1583 # As an initial value of the pixel coordinate `x_0` we take 1584 # "focal plane coordinate" `x'=W^{-1}(w)=wcs_world2pix(w)`. 1585 # We stop iterations when `|x_{i+1}-x_i|<tolerance`. We also 1586 # consider the process to be diverging if 1587 # `|x_{i+1}-x_i|>|x_i-x_{i-1}|` 1588 # **when** `|x_{i+1}-x_i|>=tolerance` (when current 1589 # approximation is close to the true solution, 1590 # `|x_{i+1}-x_i|>|x_i-x_{i-1}|` may be due to rounding errors 1591 # and we ignore such "divergences" when 1592 # `|x_{i+1}-x_i|<tolerance`). It may appear that checking for 1593 # `|x_{i+1}-x_i|<tolerance` in order to ignore divergence is 1594 # unnecessary since the iterative process should stop anyway, 1595 # however, the proposed implementation of this iterative 1596 # process is completely vectorized and, therefore, we may 1597 # continue iterating over *some* points even though they have 1598 # converged to within a specified tolerance (while iterating 1599 # over other points that have not yet converged to 1600 # a solution). 1601 # 1602 # In order to efficiently implement iterative process (6) 1603 # using available methods in `astropy.wcs.WCS`, we add and 1604 # subtract `x_i` from the right side of equation (6): 1605 # 1606 # (7) x_{i+1} = x'-(x_i+f(x_i))+x_i = x'-pix2foc(x_i)+x_i, 1607 # 1608 # where `x'=wcs_world2pix(w)` and it is computed only *once* 1609 # before the beginning of the iterative process (and we also 1610 # set `x_0=x'`). By using `pix2foc` at each iteration instead 1611 # of `all_pix2world` we get about 25% increase in performance 1612 # (by not performing the linear `W` transformation at each 1613 # step) and we also avoid the "RA wrapping" issue described 1614 # above (by working in focal plane coordinates and avoiding 1615 # pix->world transformations). 1616 # 1617 # As an added benefit, the process converges to the correct 1618 # solution in just one iteration when distortions are not 1619 # present (compare to 1620 # https://github.com/astropy/astropy/issues/1977 and 1621 # https://github.com/astropy/astropy/pull/2294): in this case 1622 # `pix2foc` is the identical transformation 1623 # `x_i=pix2foc(x_i)` and from equation (7) we get: 1624 # 1625 # x' = x_0 = wcs_world2pix(w) 1626 # x_1 = x' - pix2foc(x_0) + x_0 = x' - pix2foc(x') + x' = x' 1627 # = wcs_world2pix(w) = x_0 1628 # => 1629 # |x_1-x_0| = 0 < tolerance (with tolerance > 0) 1630 # 1631 # However, for performance reasons, it is still better to 1632 # avoid iterations altogether and return the exact linear 1633 # solution (`wcs_world2pix`) right-away when non-linear 1634 # distortions are not present by checking that attributes 1635 # `sip`, `cpdis1`, `cpdis2`, `det2im1`, and `det2im2` are 1636 # *all* `None`. 1637 # 1638 # 1639 # ### Outline of the Algorithm ### 1640 # 1641 # 1642 # While the proposed code is relatively long (considering 1643 # the simplicity of the algorithm), this is due to: 1) 1644 # checking if iterative solution is necessary at all; 2) 1645 # checking for divergence; 3) re-implementation of the 1646 # completely vectorized algorithm as an "adaptive" vectorized 1647 # algorithm (for cases when some points diverge for which we 1648 # want to stop iterations). In my tests, the adaptive version 1649 # of the algorithm is about 50% slower than non-adaptive 1650 # version for all HST images. 1651 # 1652 # The essential part of the vectorized non-adaptive algorithm 1653 # (without divergence and other checks) can be described 1654 # as follows: 1655 # 1656 # pix0 = self.wcs_world2pix(world, origin) 1657 # pix = pix0.copy() # 0-order solution 1658 # 1659 # for k in range(maxiter): 1660 # # find correction to the previous solution: 1661 # dpix = self.pix2foc(pix, origin) - pix0 1662 # 1663 # # compute norm (L2) of the correction: 1664 # dn = np.linalg.norm(dpix, axis=1) 1665 # 1666 # # apply correction: 1667 # pix -= dpix 1668 # 1669 # # check convergence: 1670 # if np.max(dn) < tolerance: 1671 # break 1672 # 1673 # return pix 1674 # 1675 # Here, the input parameter `world` can be a `MxN` array 1676 # where `M` is the number of coordinate axes in WCS and `N` 1677 # is the number of points to be converted simultaneously to 1678 # image coordinates. 1679 # 1680 # 1681 # ### IMPORTANT NOTE: ### 1682 # 1683 # If, in the future releases of the `~astropy.wcs`, 1684 # `pix2foc` will not apply all the required distortion 1685 # corrections then in the code below, calls to `pix2foc` will 1686 # have to be replaced with 1687 # wcs_world2pix(all_pix2world(pix_list, origin), origin) 1688 # 1689 1690 # ############################################################ 1691 # # INITIALIZE ITERATIVE PROCESS: ## 1692 # ############################################################ 1693 1694 # initial approximation (linear WCS based only) 1695 pix0 = self.wcs_world2pix(world, origin) 1696 1697 # Check that an iterative solution is required at all 1698 # (when any of the non-CD-matrix-based corrections are 1699 # present). If not required return the initial 1700 # approximation (pix0). 1701 if not self.has_distortion: 1702 # No non-WCS corrections detected so 1703 # simply return initial approximation: 1704 return pix0 1705 1706 pix = pix0.copy() # 0-order solution 1707 1708 # initial correction: 1709 dpix = self.pix2foc(pix, origin) - pix0 1710 1711 # Update initial solution: 1712 pix -= dpix 1713 1714 # Norm (L2) squared of the correction: 1715 dn = np.sum(dpix*dpix, axis=1) 1716 dnprev = dn.copy() # if adaptive else dn 1717 tol2 = tolerance**2 1718 1719 # Prepare for iterative process 1720 k = 1 1721 ind = None 1722 inddiv = None 1723 1724 # Turn off numpy runtime warnings for 'invalid' and 'over': 1725 old_invalid = np.geterr()['invalid'] 1726 old_over = np.geterr()['over'] 1727 np.seterr(invalid='ignore', over='ignore') 1728 1729 # ############################################################ 1730 # # NON-ADAPTIVE ITERATIONS: ## 1731 # ############################################################ 1732 if not adaptive: 1733 # Fixed-point iterations: 1734 while (np.nanmax(dn) >= tol2 and k < maxiter): 1735 # Find correction to the previous solution: 1736 dpix = self.pix2foc(pix, origin) - pix0 1737 1738 # Compute norm (L2) squared of the correction: 1739 dn = np.sum(dpix*dpix, axis=1) 1740 1741 # Check for divergence (we do this in two stages 1742 # to optimize performance for the most common 1743 # scenario when successive approximations converge): 1744 if detect_divergence: 1745 divergent = (dn >= dnprev) 1746 if np.any(divergent): 1747 # Find solutions that have not yet converged: 1748 slowconv = (dn >= tol2) 1749 inddiv, = np.where(divergent & slowconv) 1750 1751 if inddiv.shape[0] > 0: 1752 # Update indices of elements that 1753 # still need correction: 1754 conv = (dn < dnprev) 1755 iconv = np.where(conv) 1756 1757 # Apply correction: 1758 dpixgood = dpix[iconv] 1759 pix[iconv] -= dpixgood 1760 dpix[iconv] = dpixgood 1761 1762 # For the next iteration choose 1763 # non-divergent points that have not yet 1764 # converged to the requested accuracy: 1765 ind, = np.where(slowconv & conv) 1766 pix0 = pix0[ind] 1767 dnprev[ind] = dn[ind] 1768 k += 1 1769 1770 # Switch to adaptive iterations: 1771 adaptive = True 1772 break 1773 # Save current correction magnitudes for later: 1774 dnprev = dn 1775 1776 # Apply correction: 1777 pix -= dpix 1778 k += 1 1779 1780 # ############################################################ 1781 # # ADAPTIVE ITERATIONS: ## 1782 # ############################################################ 1783 if adaptive: 1784 if ind is None: 1785 ind, = np.where(np.isfinite(pix).all(axis=1)) 1786 pix0 = pix0[ind] 1787 1788 # "Adaptive" fixed-point iterations: 1789 while (ind.shape[0] > 0 and k < maxiter): 1790 # Find correction to the previous solution: 1791 dpixnew = self.pix2foc(pix[ind], origin) - pix0 1792 1793 # Compute norm (L2) of the correction: 1794 dnnew = np.sum(np.square(dpixnew), axis=1) 1795 1796 # Bookkeeping of corrections: 1797 dnprev[ind] = dn[ind].copy() 1798 dn[ind] = dnnew 1799 1800 if detect_divergence: 1801 # Find indices of pixels that are converging: 1802 conv = (dnnew < dnprev[ind]) 1803 iconv = np.where(conv) 1804 iiconv = ind[iconv] 1805 1806 # Apply correction: 1807 dpixgood = dpixnew[iconv] 1808 pix[iiconv] -= dpixgood 1809 dpix[iiconv] = dpixgood 1810 1811 # Find indices of solutions that have not yet 1812 # converged to the requested accuracy 1813 # AND that do not diverge: 1814 subind, = np.where((dnnew >= tol2) & conv) 1815 1816 else: 1817 # Apply correction: 1818 pix[ind] -= dpixnew 1819 dpix[ind] = dpixnew 1820 1821 # Find indices of solutions that have not yet 1822 # converged to the requested accuracy: 1823 subind, = np.where(dnnew >= tol2) 1824 1825 # Choose solutions that need more iterations: 1826 ind = ind[subind] 1827 pix0 = pix0[subind] 1828 1829 k += 1 1830 1831 # ############################################################ 1832 # # FINAL DETECTION OF INVALID, DIVERGING, ## 1833 # # AND FAILED-TO-CONVERGE POINTS ## 1834 # ############################################################ 1835 # Identify diverging and/or invalid points: 1836 invalid = ((~np.all(np.isfinite(pix), axis=1)) & 1837 (np.all(np.isfinite(world), axis=1))) 1838 1839 # When detect_divergence==False, dnprev is outdated 1840 # (it is the norm of the very first correction). 1841 # Still better than nothing... 1842 inddiv, = np.where(((dn >= tol2) & (dn >= dnprev)) | invalid) 1843 if inddiv.shape[0] == 0: 1844 inddiv = None 1845 1846 # Identify points that did not converge within 'maxiter' 1847 # iterations: 1848 if k >= maxiter: 1849 ind, = np.where((dn >= tol2) & (dn < dnprev) & (~invalid)) 1850 if ind.shape[0] == 0: 1851 ind = None 1852 else: 1853 ind = None 1854 1855 # Restore previous numpy error settings: 1856 np.seterr(invalid=old_invalid, over=old_over) 1857 1858 # ############################################################ 1859 # # RAISE EXCEPTION IF DIVERGING OR TOO SLOWLY CONVERGING ## 1860 # # DATA POINTS HAVE BEEN DETECTED: ## 1861 # ############################################################ 1862 if (ind is not None or inddiv is not None) and not quiet: 1863 if inddiv is None: 1864 raise NoConvergence( 1865 "'WCS.all_world2pix' failed to " 1866 "converge to the requested accuracy after {:d} " 1867 "iterations.".format(k), best_solution=pix, 1868 accuracy=np.abs(dpix), niter=k, 1869 slow_conv=ind, divergent=None) 1870 else: 1871 raise NoConvergence( 1872 "'WCS.all_world2pix' failed to " 1873 "converge to the requested accuracy.\n" 1874 "After {:d} iterations, the solution is diverging " 1875 "at least for one input point." 1876 .format(k), best_solution=pix, 1877 accuracy=np.abs(dpix), niter=k, 1878 slow_conv=ind, divergent=inddiv) 1879 1880 return pix 1881 1882 @deprecated_renamed_argument('accuracy', 'tolerance', '4.3') 1883 def all_world2pix(self, *args, tolerance=1e-4, maxiter=20, adaptive=False, 1884 detect_divergence=True, quiet=False, **kwargs): 1885 if self.wcs is None: 1886 raise ValueError("No basic WCS settings were created.") 1887 1888 return self._array_converter( 1889 lambda *args, **kwargs: 1890 self._all_world2pix( 1891 *args, tolerance=tolerance, maxiter=maxiter, 1892 adaptive=adaptive, detect_divergence=detect_divergence, 1893 quiet=quiet), 1894 'input', *args, **kwargs 1895 ) 1896 1897 all_world2pix.__doc__ = """ 1898 all_world2pix(*arg, tolerance=1.0e-4, maxiter=20, 1899 adaptive=False, detect_divergence=True, quiet=False) 1900 1901 Transforms world coordinates to pixel coordinates, using 1902 numerical iteration to invert the full forward transformation 1903 `~astropy.wcs.WCS.all_pix2world` with complete 1904 distortion model. 1905 1906 1907 Parameters 1908 ---------- 1909 {0} 1910 1911 For a transformation that is not two-dimensional, the 1912 two-argument form must be used. 1913 1914 {1} 1915 1916 tolerance : float, optional (default = 1.0e-4) 1917 Tolerance of solution. Iteration terminates when the 1918 iterative solver estimates that the "true solution" is 1919 within this many pixels current estimate, more 1920 specifically, when the correction to the solution found 1921 during the previous iteration is smaller 1922 (in the sense of the L2 norm) than ``tolerance``. 1923 1924 maxiter : int, optional (default = 20) 1925 Maximum number of iterations allowed to reach a solution. 1926 1927 quiet : bool, optional (default = False) 1928 Do not throw :py:class:`NoConvergence` exceptions when 1929 the method does not converge to a solution with the 1930 required accuracy within a specified number of maximum 1931 iterations set by ``maxiter`` parameter. Instead, 1932 simply return the found solution. 1933 1934 Other Parameters 1935 ---------------- 1936 adaptive : bool, optional (default = False) 1937 Specifies whether to adaptively select only points that 1938 did not converge to a solution within the required 1939 accuracy for the next iteration. Default is recommended 1940 for HST as well as most other instruments. 1941 1942 .. note:: 1943 The :py:meth:`all_world2pix` uses a vectorized 1944 implementation of the method of consecutive 1945 approximations (see ``Notes`` section below) in which it 1946 iterates over *all* input points *regardless* until 1947 the required accuracy has been reached for *all* input 1948 points. In some cases it may be possible that 1949 *almost all* points have reached the required accuracy 1950 but there are only a few of input data points for 1951 which additional iterations may be needed (this 1952 depends mostly on the characteristics of the geometric 1953 distortions for a given instrument). In this situation 1954 it may be advantageous to set ``adaptive`` = `True` in 1955 which case :py:meth:`all_world2pix` will continue 1956 iterating *only* over the points that have not yet 1957 converged to the required accuracy. However, for the 1958 HST's ACS/WFC detector, which has the strongest 1959 distortions of all HST instruments, testing has 1960 shown that enabling this option would lead to a about 1961 50-100% penalty in computational time (depending on 1962 specifics of the image, geometric distortions, and 1963 number of input points to be converted). Therefore, 1964 for HST and possibly instruments, it is recommended 1965 to set ``adaptive`` = `False`. The only danger in 1966 getting this setting wrong will be a performance 1967 penalty. 1968 1969 .. note:: 1970 When ``detect_divergence`` is `True`, 1971 :py:meth:`all_world2pix` will automatically switch 1972 to the adaptive algorithm once divergence has been 1973 detected. 1974 1975 detect_divergence : bool, optional (default = True) 1976 Specifies whether to perform a more detailed analysis 1977 of the convergence to a solution. Normally 1978 :py:meth:`all_world2pix` may not achieve the required 1979 accuracy if either the ``tolerance`` or ``maxiter`` arguments 1980 are too low. However, it may happen that for some 1981 geometric distortions the conditions of convergence for 1982 the the method of consecutive approximations used by 1983 :py:meth:`all_world2pix` may not be satisfied, in which 1984 case consecutive approximations to the solution will 1985 diverge regardless of the ``tolerance`` or ``maxiter`` 1986 settings. 1987 1988 When ``detect_divergence`` is `False`, these divergent 1989 points will be detected as not having achieved the 1990 required accuracy (without further details). In addition, 1991 if ``adaptive`` is `False` then the algorithm will not 1992 know that the solution (for specific points) is diverging 1993 and will continue iterating and trying to "improve" 1994 diverging solutions. This may result in ``NaN`` or 1995 ``Inf`` values in the return results (in addition to a 1996 performance penalties). Even when ``detect_divergence`` 1997 is `False`, :py:meth:`all_world2pix`, at the end of the 1998 iterative process, will identify invalid results 1999 (``NaN`` or ``Inf``) as "diverging" solutions and will 2000 raise :py:class:`NoConvergence` unless the ``quiet`` 2001 parameter is set to `True`. 2002 2003 When ``detect_divergence`` is `True`, 2004 :py:meth:`all_world2pix` will detect points for which 2005 current correction to the coordinates is larger than 2006 the correction applied during the previous iteration 2007 **if** the requested accuracy **has not yet been 2008 achieved**. In this case, if ``adaptive`` is `True`, 2009 these points will be excluded from further iterations and 2010 if ``adaptive`` is `False`, :py:meth:`all_world2pix` will 2011 automatically switch to the adaptive algorithm. Thus, the 2012 reported divergent solution will be the latest converging 2013 solution computed immediately *before* divergence 2014 has been detected. 2015 2016 .. note:: 2017 When accuracy has been achieved, small increases in 2018 current corrections may be possible due to rounding 2019 errors (when ``adaptive`` is `False`) and such 2020 increases will be ignored. 2021 2022 .. note:: 2023 Based on our testing using HST ACS/WFC images, setting 2024 ``detect_divergence`` to `True` will incur about 5-20% 2025 performance penalty with the larger penalty 2026 corresponding to ``adaptive`` set to `True`. 2027 Because the benefits of enabling this 2028 feature outweigh the small performance penalty, 2029 especially when ``adaptive`` = `False`, it is 2030 recommended to set ``detect_divergence`` to `True`, 2031 unless extensive testing of the distortion models for 2032 images from specific instruments show a good stability 2033 of the numerical method for a wide range of 2034 coordinates (even outside the image itself). 2035 2036 .. note:: 2037 Indices of the diverging inverse solutions will be 2038 reported in the ``divergent`` attribute of the 2039 raised :py:class:`NoConvergence` exception object. 2040 2041 Returns 2042 ------- 2043 2044 {2} 2045 2046 Notes 2047 ----- 2048 The order of the axes for the input world array is determined by 2049 the ``CTYPEia`` keywords in the FITS header, therefore it may 2050 not always be of the form (*ra*, *dec*). The 2051 `~astropy.wcs.Wcsprm.lat`, `~astropy.wcs.Wcsprm.lng`, 2052 `~astropy.wcs.Wcsprm.lattyp`, and 2053 `~astropy.wcs.Wcsprm.lngtyp` 2054 members can be used to determine the order of the axes. 2055 2056 Using the method of fixed-point iterations approximations we 2057 iterate starting with the initial approximation, which is 2058 computed using the non-distortion-aware 2059 :py:meth:`wcs_world2pix` (or equivalent). 2060 2061 The :py:meth:`all_world2pix` function uses a vectorized 2062 implementation of the method of consecutive approximations and 2063 therefore it is highly efficient (>30x) when *all* data points 2064 that need to be converted from sky coordinates to image 2065 coordinates are passed at *once*. Therefore, it is advisable, 2066 whenever possible, to pass as input a long array of all points 2067 that need to be converted to :py:meth:`all_world2pix` instead 2068 of calling :py:meth:`all_world2pix` for each data point. Also 2069 see the note to the ``adaptive`` parameter. 2070 2071 Raises 2072 ------ 2073 NoConvergence 2074 The method did not converge to a 2075 solution to the required accuracy within a specified 2076 number of maximum iterations set by the ``maxiter`` 2077 parameter. To turn off this exception, set ``quiet`` to 2078 `True`. Indices of the points for which the requested 2079 accuracy was not achieved (if any) will be listed in the 2080 ``slow_conv`` attribute of the 2081 raised :py:class:`NoConvergence` exception object. 2082 2083 See :py:class:`NoConvergence` documentation for 2084 more details. 2085 2086 MemoryError 2087 Memory allocation failed. 2088 2089 SingularMatrixError 2090 Linear transformation matrix is singular. 2091 2092 InconsistentAxisTypesError 2093 Inconsistent or unrecognized coordinate axis types. 2094 2095 ValueError 2096 Invalid parameter value. 2097 2098 ValueError 2099 Invalid coordinate transformation parameters. 2100 2101 ValueError 2102 x- and y-coordinate arrays are not the same size. 2103 2104 InvalidTransformError 2105 Invalid coordinate transformation parameters. 2106 2107 InvalidTransformError 2108 Ill-conditioned coordinate transformation parameters. 2109 2110 Examples 2111 -------- 2112 >>> import astropy.io.fits as fits 2113 >>> import astropy.wcs as wcs 2114 >>> import numpy as np 2115 >>> import os 2116 2117 >>> filename = os.path.join(wcs.__path__[0], 'tests/data/j94f05bgq_flt.fits') 2118 >>> hdulist = fits.open(filename) 2119 >>> w = wcs.WCS(hdulist[('sci',1)].header, hdulist) 2120 >>> hdulist.close() 2121 2122 >>> ra, dec = w.all_pix2world([1,2,3], [1,1,1], 1) 2123 >>> print(ra) # doctest: +FLOAT_CMP 2124 [ 5.52645627 5.52649663 5.52653698] 2125 >>> print(dec) # doctest: +FLOAT_CMP 2126 [-72.05171757 -72.05171276 -72.05170795] 2127 >>> radec = w.all_pix2world([[1,1], [2,1], [3,1]], 1) 2128 >>> print(radec) # doctest: +FLOAT_CMP 2129 [[ 5.52645627 -72.05171757] 2130 [ 5.52649663 -72.05171276] 2131 [ 5.52653698 -72.05170795]] 2132 >>> x, y = w.all_world2pix(ra, dec, 1) 2133 >>> print(x) # doctest: +FLOAT_CMP 2134 [ 1.00000238 2.00000237 3.00000236] 2135 >>> print(y) # doctest: +FLOAT_CMP 2136 [ 0.99999996 0.99999997 0.99999997] 2137 >>> xy = w.all_world2pix(radec, 1) 2138 >>> print(xy) # doctest: +FLOAT_CMP 2139 [[ 1.00000238 0.99999996] 2140 [ 2.00000237 0.99999997] 2141 [ 3.00000236 0.99999997]] 2142 >>> xy = w.all_world2pix(radec, 1, maxiter=3, 2143 ... tolerance=1.0e-10, quiet=False) 2144 Traceback (most recent call last): 2145 ... 2146 NoConvergence: 'WCS.all_world2pix' failed to converge to the 2147 requested accuracy. After 3 iterations, the solution is 2148 diverging at least for one input point. 2149 2150 >>> # Now try to use some diverging data: 2151 >>> divradec = w.all_pix2world([[1.0, 1.0], 2152 ... [10000.0, 50000.0], 2153 ... [3.0, 1.0]], 1) 2154 >>> print(divradec) # doctest: +FLOAT_CMP 2155 [[ 5.52645627 -72.05171757] 2156 [ 7.15976932 -70.8140779 ] 2157 [ 5.52653698 -72.05170795]] 2158 2159 >>> # First, turn detect_divergence on: 2160 >>> try: # doctest: +FLOAT_CMP 2161 ... xy = w.all_world2pix(divradec, 1, maxiter=20, 2162 ... tolerance=1.0e-4, adaptive=False, 2163 ... detect_divergence=True, 2164 ... quiet=False) 2165 ... except wcs.wcs.NoConvergence as e: 2166 ... print("Indices of diverging points: {{0}}" 2167 ... .format(e.divergent)) 2168 ... print("Indices of poorly converging points: {{0}}" 2169 ... .format(e.slow_conv)) 2170 ... print("Best solution:\\n{{0}}".format(e.best_solution)) 2171 ... print("Achieved accuracy:\\n{{0}}".format(e.accuracy)) 2172 Indices of diverging points: [1] 2173 Indices of poorly converging points: None 2174 Best solution: 2175 [[ 1.00000238e+00 9.99999965e-01] 2176 [ -1.99441636e+06 1.44309097e+06] 2177 [ 3.00000236e+00 9.99999966e-01]] 2178 Achieved accuracy: 2179 [[ 6.13968380e-05 8.59638593e-07] 2180 [ 8.59526812e+11 6.61713548e+11] 2181 [ 6.09398446e-05 8.38759724e-07]] 2182 >>> raise e 2183 Traceback (most recent call last): 2184 ... 2185 NoConvergence: 'WCS.all_world2pix' failed to converge to the 2186 requested accuracy. After 5 iterations, the solution is 2187 diverging at least for one input point. 2188 2189 >>> # This time turn detect_divergence off: 2190 >>> try: # doctest: +FLOAT_CMP 2191 ... xy = w.all_world2pix(divradec, 1, maxiter=20, 2192 ... tolerance=1.0e-4, adaptive=False, 2193 ... detect_divergence=False, 2194 ... quiet=False) 2195 ... except wcs.wcs.NoConvergence as e: 2196 ... print("Indices of diverging points: {{0}}" 2197 ... .format(e.divergent)) 2198 ... print("Indices of poorly converging points: {{0}}" 2199 ... .format(e.slow_conv)) 2200 ... print("Best solution:\\n{{0}}".format(e.best_solution)) 2201 ... print("Achieved accuracy:\\n{{0}}".format(e.accuracy)) 2202 Indices of diverging points: [1] 2203 Indices of poorly converging points: None 2204 Best solution: 2205 [[ 1.00000009 1. ] 2206 [ nan nan] 2207 [ 3.00000009 1. ]] 2208 Achieved accuracy: 2209 [[ 2.29417358e-06 3.21222995e-08] 2210 [ nan nan] 2211 [ 2.27407877e-06 3.13005639e-08]] 2212 >>> raise e 2213 Traceback (most recent call last): 2214 ... 2215 NoConvergence: 'WCS.all_world2pix' failed to converge to the 2216 requested accuracy. After 6 iterations, the solution is 2217 diverging at least for one input point. 2218 2219 """.format(docstrings.TWO_OR_MORE_ARGS('naxis', 8), 2220 docstrings.RA_DEC_ORDER(8), 2221 docstrings.RETURNS('pixel coordinates', 8)) 2222 2223 def wcs_world2pix(self, *args, **kwargs): 2224 if self.wcs is None: 2225 raise ValueError("No basic WCS settings were created.") 2226 return self._array_converter( 2227 lambda xy, o: self.wcs.s2p(xy, o)['pixcrd'], 2228 'input', *args, **kwargs) 2229 wcs_world2pix.__doc__ = """ 2230 Transforms world coordinates to pixel coordinates, using only 2231 the basic `wcslib`_ WCS transformation. No `SIP`_ or 2232 `distortion paper`_ table lookup transformation is applied. 2233 2234 Parameters 2235 ---------- 2236 {} 2237 2238 For a transformation that is not two-dimensional, the 2239 two-argument form must be used. 2240 2241 {} 2242 2243 Returns 2244 ------- 2245 2246 {} 2247 2248 Notes 2249 ----- 2250 The order of the axes for the input world array is determined by 2251 the ``CTYPEia`` keywords in the FITS header, therefore it may 2252 not always be of the form (*ra*, *dec*). The 2253 `~astropy.wcs.Wcsprm.lat`, `~astropy.wcs.Wcsprm.lng`, 2254 `~astropy.wcs.Wcsprm.lattyp` and `~astropy.wcs.Wcsprm.lngtyp` 2255 members can be used to determine the order of the axes. 2256 2257 Raises 2258 ------ 2259 MemoryError 2260 Memory allocation failed. 2261 2262 SingularMatrixError 2263 Linear transformation matrix is singular. 2264 2265 InconsistentAxisTypesError 2266 Inconsistent or unrecognized coordinate axis types. 2267 2268 ValueError 2269 Invalid parameter value. 2270 2271 ValueError 2272 Invalid coordinate transformation parameters. 2273 2274 ValueError 2275 x- and y-coordinate arrays are not the same size. 2276 2277 InvalidTransformError 2278 Invalid coordinate transformation parameters. 2279 2280 InvalidTransformError 2281 Ill-conditioned coordinate transformation parameters. 2282 """.format(docstrings.TWO_OR_MORE_ARGS('naxis', 8), 2283 docstrings.RA_DEC_ORDER(8), 2284 docstrings.RETURNS('pixel coordinates', 8)) 2285 2286 def pix2foc(self, *args): 2287 return self._array_converter(self._pix2foc, None, *args) 2288 pix2foc.__doc__ = """ 2289 Convert pixel coordinates to focal plane coordinates using the 2290 `SIP`_ polynomial distortion convention and `distortion 2291 paper`_ table-lookup correction. 2292 2293 The output is in absolute pixel coordinates, not relative to 2294 ``CRPIX``. 2295 2296 Parameters 2297 ---------- 2298 2299 {} 2300 2301 Returns 2302 ------- 2303 2304 {} 2305 2306 Raises 2307 ------ 2308 MemoryError 2309 Memory allocation failed. 2310 2311 ValueError 2312 Invalid coordinate transformation parameters. 2313 """.format(docstrings.TWO_OR_MORE_ARGS('2', 8), 2314 docstrings.RETURNS('focal coordinates', 8)) 2315 2316 def p4_pix2foc(self, *args): 2317 return self._array_converter(self._p4_pix2foc, None, *args) 2318 p4_pix2foc.__doc__ = """ 2319 Convert pixel coordinates to focal plane coordinates using 2320 `distortion paper`_ table-lookup correction. 2321 2322 The output is in absolute pixel coordinates, not relative to 2323 ``CRPIX``. 2324 2325 Parameters 2326 ---------- 2327 2328 {} 2329 2330 Returns 2331 ------- 2332 2333 {} 2334 2335 Raises 2336 ------ 2337 MemoryError 2338 Memory allocation failed. 2339 2340 ValueError 2341 Invalid coordinate transformation parameters. 2342 """.format(docstrings.TWO_OR_MORE_ARGS('2', 8), 2343 docstrings.RETURNS('focal coordinates', 8)) 2344 2345 def det2im(self, *args): 2346 return self._array_converter(self._det2im, None, *args) 2347 det2im.__doc__ = """ 2348 Convert detector coordinates to image plane coordinates using 2349 `distortion paper`_ table-lookup correction. 2350 2351 The output is in absolute pixel coordinates, not relative to 2352 ``CRPIX``. 2353 2354 Parameters 2355 ---------- 2356 2357 {} 2358 2359 Returns 2360 ------- 2361 2362 {} 2363 2364 Raises 2365 ------ 2366 MemoryError 2367 Memory allocation failed. 2368 2369 ValueError 2370 Invalid coordinate transformation parameters. 2371 """.format(docstrings.TWO_OR_MORE_ARGS('2', 8), 2372 docstrings.RETURNS('pixel coordinates', 8)) 2373 2374 def sip_pix2foc(self, *args): 2375 if self.sip is None: 2376 if len(args) == 2: 2377 return args[0] 2378 elif len(args) == 3: 2379 return args[:2] 2380 else: 2381 raise TypeError("Wrong number of arguments") 2382 return self._array_converter(self.sip.pix2foc, None, *args) 2383 sip_pix2foc.__doc__ = """ 2384 Convert pixel coordinates to focal plane coordinates using the 2385 `SIP`_ polynomial distortion convention. 2386 2387 The output is in pixel coordinates, relative to ``CRPIX``. 2388 2389 FITS WCS `distortion paper`_ table lookup correction is not 2390 applied, even if that information existed in the FITS file 2391 that initialized this :class:`~astropy.wcs.WCS` object. To 2392 correct for that, use `~astropy.wcs.WCS.pix2foc` or 2393 `~astropy.wcs.WCS.p4_pix2foc`. 2394 2395 Parameters 2396 ---------- 2397 2398 {} 2399 2400 Returns 2401 ------- 2402 2403 {} 2404 2405 Raises 2406 ------ 2407 MemoryError 2408 Memory allocation failed. 2409 2410 ValueError 2411 Invalid coordinate transformation parameters. 2412 """.format(docstrings.TWO_OR_MORE_ARGS('2', 8), 2413 docstrings.RETURNS('focal coordinates', 8)) 2414 2415 def sip_foc2pix(self, *args): 2416 if self.sip is None: 2417 if len(args) == 2: 2418 return args[0] 2419 elif len(args) == 3: 2420 return args[:2] 2421 else: 2422 raise TypeError("Wrong number of arguments") 2423 return self._array_converter(self.sip.foc2pix, None, *args) 2424 sip_foc2pix.__doc__ = """ 2425 Convert focal plane coordinates to pixel coordinates using the 2426 `SIP`_ polynomial distortion convention. 2427 2428 FITS WCS `distortion paper`_ table lookup distortion 2429 correction is not applied, even if that information existed in 2430 the FITS file that initialized this `~astropy.wcs.WCS` object. 2431 2432 Parameters 2433 ---------- 2434 2435 {} 2436 2437 Returns 2438 ------- 2439 2440 {} 2441 2442 Raises 2443 ------ 2444 MemoryError 2445 Memory allocation failed. 2446 2447 ValueError 2448 Invalid coordinate transformation parameters. 2449 """.format(docstrings.TWO_OR_MORE_ARGS('2', 8), 2450 docstrings.RETURNS('pixel coordinates', 8)) 2451 2452 def proj_plane_pixel_scales(self): 2453 """ 2454 Calculate pixel scales along each axis of the image pixel at 2455 the ``CRPIX`` location once it is projected onto the 2456 "plane of intermediate world coordinates" as defined in 2457 `Greisen & Calabretta 2002, A&A, 395, 1061 <https://ui.adsabs.harvard.edu/abs/2002A%26A...395.1061G>`_. 2458 2459 .. note:: 2460 This method is concerned **only** about the transformation 2461 "image plane"->"projection plane" and **not** about the 2462 transformation "celestial sphere"->"projection plane"->"image plane". 2463 Therefore, this function ignores distortions arising due to 2464 non-linear nature of most projections. 2465 2466 .. note:: 2467 This method only returns sensible answers if the WCS contains 2468 celestial axes, i.e., the `~astropy.wcs.WCS.celestial` WCS object. 2469 2470 Returns 2471 ------- 2472 scale : list of `~astropy.units.Quantity` 2473 A vector of projection plane increments corresponding to each 2474 pixel side (axis). 2475 2476 See Also 2477 -------- 2478 astropy.wcs.utils.proj_plane_pixel_scales 2479 2480 """ # noqa: E501 2481 from astropy.wcs.utils import proj_plane_pixel_scales # Avoid circular import 2482 values = proj_plane_pixel_scales(self) 2483 units = [u.Unit(x) for x in self.wcs.cunit] 2484 return [value * unit for (value, unit) in zip(values, units)] # Can have different units 2485 2486 def proj_plane_pixel_area(self): 2487 """ 2488 For a **celestial** WCS (see `astropy.wcs.WCS.celestial`), returns pixel 2489 area of the image pixel at the ``CRPIX`` location once it is projected 2490 onto the "plane of intermediate world coordinates" as defined in 2491 `Greisen & Calabretta 2002, A&A, 395, 1061 <https://ui.adsabs.harvard.edu/abs/2002A%26A...395.1061G>`_. 2492 2493 .. note:: 2494 This function is concerned **only** about the transformation 2495 "image plane"->"projection plane" and **not** about the 2496 transformation "celestial sphere"->"projection plane"->"image plane". 2497 Therefore, this function ignores distortions arising due to 2498 non-linear nature of most projections. 2499 2500 .. note:: 2501 This method only returns sensible answers if the WCS contains 2502 celestial axes, i.e., the `~astropy.wcs.WCS.celestial` WCS object. 2503 2504 Returns 2505 ------- 2506 area : `~astropy.units.Quantity` 2507 Area (in the projection plane) of the pixel at ``CRPIX`` location. 2508 2509 Raises 2510 ------ 2511 ValueError 2512 Pixel area is defined only for 2D pixels. Most likely the 2513 `~astropy.wcs.Wcsprm.cd` matrix of the `~astropy.wcs.WCS.celestial` 2514 WCS is not a square matrix of second order. 2515 2516 Notes 2517 ----- 2518 2519 Depending on the application, square root of the pixel area can be used to 2520 represent a single pixel scale of an equivalent square pixel 2521 whose area is equal to the area of a generally non-square pixel. 2522 2523 See Also 2524 -------- 2525 astropy.wcs.utils.proj_plane_pixel_area 2526 2527 """ # noqa: E501 2528 from astropy.wcs.utils import proj_plane_pixel_area # Avoid circular import 2529 value = proj_plane_pixel_area(self) 2530 unit = u.Unit(self.wcs.cunit[0]) * u.Unit(self.wcs.cunit[1]) # 2D only 2531 return value * unit 2532 2533 def to_fits(self, relax=False, key=None): 2534 """ 2535 Generate an `~astropy.io.fits.HDUList` object with all of the 2536 information stored in this object. This should be logically identical 2537 to the input FITS file, but it will be normalized in a number of ways. 2538 2539 See `to_header` for some warnings about the output produced. 2540 2541 Parameters 2542 ---------- 2543 2544 relax : bool or int, optional 2545 Degree of permissiveness: 2546 2547 - `False` (default): Write all extensions that are 2548 considered to be safe and recommended. 2549 2550 - `True`: Write all recognized informal extensions of the 2551 WCS standard. 2552 2553 - `int`: a bit field selecting specific extensions to 2554 write. See :ref:`astropy:relaxwrite` for details. 2555 2556 key : str 2557 The name of a particular WCS transform to use. This may be 2558 either ``' '`` or ``'A'``-``'Z'`` and corresponds to the ``"a"`` 2559 part of the ``CTYPEia`` cards. 2560 2561 Returns 2562 ------- 2563 hdulist : `~astropy.io.fits.HDUList` 2564 """ 2565 2566 header = self.to_header(relax=relax, key=key) 2567 2568 hdu = fits.PrimaryHDU(header=header) 2569 hdulist = fits.HDUList(hdu) 2570 2571 self._write_det2im(hdulist) 2572 self._write_distortion_kw(hdulist) 2573 2574 return hdulist 2575 2576 def to_header(self, relax=None, key=None): 2577 """Generate an `astropy.io.fits.Header` object with the basic WCS 2578 and SIP information stored in this object. This should be 2579 logically identical to the input FITS file, but it will be 2580 normalized in a number of ways. 2581 2582 .. warning:: 2583 2584 This function does not write out FITS WCS `distortion 2585 paper`_ information, since that requires multiple FITS 2586 header data units. To get a full representation of 2587 everything in this object, use `to_fits`. 2588 2589 Parameters 2590 ---------- 2591 relax : bool or int, optional 2592 Degree of permissiveness: 2593 2594 - `False` (default): Write all extensions that are 2595 considered to be safe and recommended. 2596 2597 - `True`: Write all recognized informal extensions of the 2598 WCS standard. 2599 2600 - `int`: a bit field selecting specific extensions to 2601 write. See :ref:`astropy:relaxwrite` for details. 2602 2603 If the ``relax`` keyword argument is not given and any 2604 keywords were omitted from the output, an 2605 `~astropy.utils.exceptions.AstropyWarning` is displayed. 2606 To override this, explicitly pass a value to ``relax``. 2607 2608 key : str 2609 The name of a particular WCS transform to use. This may be 2610 either ``' '`` or ``'A'``-``'Z'`` and corresponds to the ``"a"`` 2611 part of the ``CTYPEia`` cards. 2612 2613 Returns 2614 ------- 2615 header : `astropy.io.fits.Header` 2616 2617 Notes 2618 ----- 2619 The output header will almost certainly differ from the input in a 2620 number of respects: 2621 2622 1. The output header only contains WCS-related keywords. In 2623 particular, it does not contain syntactically-required 2624 keywords such as ``SIMPLE``, ``NAXIS``, ``BITPIX``, or 2625 ``END``. 2626 2627 2. Deprecated (e.g. ``CROTAn``) or non-standard usage will 2628 be translated to standard (this is partially dependent on 2629 whether ``fix`` was applied). 2630 2631 3. Quantities will be converted to the units used internally, 2632 basically SI with the addition of degrees. 2633 2634 4. Floating-point quantities may be given to a different decimal 2635 precision. 2636 2637 5. Elements of the ``PCi_j`` matrix will be written if and 2638 only if they differ from the unit matrix. Thus, if the 2639 matrix is unity then no elements will be written. 2640 2641 6. Additional keywords such as ``WCSAXES``, ``CUNITia``, 2642 ``LONPOLEa`` and ``LATPOLEa`` may appear. 2643 2644 7. The original keycomments will be lost, although 2645 `to_header` tries hard to write meaningful comments. 2646 2647 8. Keyword order may be changed. 2648 2649 """ 2650 # default precision for numerical WCS keywords 2651 precision = WCSHDO_P14 # Defined by C-ext # noqa: F821 2652 display_warning = False 2653 if relax is None: 2654 display_warning = True 2655 relax = False 2656 2657 if relax not in (True, False): 2658 do_sip = relax & WCSHDO_SIP 2659 relax &= ~WCSHDO_SIP 2660 else: 2661 do_sip = relax 2662 relax = WCSHDO_all if relax is True else WCSHDO_safe # Defined by C-ext # noqa: F821 2663 2664 relax = precision | relax 2665 2666 if self.wcs is not None: 2667 if key is not None: 2668 orig_key = self.wcs.alt 2669 self.wcs.alt = key 2670 header_string = self.wcs.to_header(relax) 2671 header = fits.Header.fromstring(header_string) 2672 keys_to_remove = ["", " ", "COMMENT"] 2673 for kw in keys_to_remove: 2674 if kw in header: 2675 del header[kw] 2676 # Check if we can handle TPD distortion correctly 2677 if int(_parsed_version[0]) * 10 + int(_parsed_version[1]) < 71: 2678 for kw, val in header.items(): 2679 if kw[:5] in ('CPDIS', 'CQDIS') and val == 'TPD': 2680 warnings.warn( 2681 f"WCS contains a TPD distortion model in {kw}. WCSLIB " 2682 f"{_wcs.__version__} is writing this in a format incompatible with " 2683 f"current versions - please update to 7.4 or use the bundled WCSLIB.", 2684 AstropyWarning) 2685 elif int(_parsed_version[0]) * 10 + int(_parsed_version[1]) < 74: 2686 for kw, val in header.items(): 2687 if kw[:5] in ('CPDIS', 'CQDIS') and val == 'TPD': 2688 warnings.warn( 2689 f"WCS contains a TPD distortion model in {kw}, which requires WCSLIB " 2690 f"7.4 or later to store in a FITS header (having {_wcs.__version__}).", 2691 AstropyWarning) 2692 else: 2693 header = fits.Header() 2694 2695 if do_sip and self.sip is not None: 2696 if self.wcs is not None and any(not ctyp.endswith('-SIP') for ctyp in self.wcs.ctype): 2697 self._fix_ctype(header, add_sip=True) 2698 2699 for kw, val in self._write_sip_kw().items(): 2700 header[kw] = val 2701 2702 if not do_sip and self.wcs is not None and any(self.wcs.ctype) and self.sip is not None: 2703 # This is called when relax is not False or WCSHDO_SIP 2704 # The default case of ``relax=None`` is handled further in the code. 2705 header = self._fix_ctype(header, add_sip=False) 2706 2707 if display_warning: 2708 full_header = self.to_header(relax=True, key=key) 2709 missing_keys = [] 2710 for kw, val in full_header.items(): 2711 if kw not in header: 2712 missing_keys.append(kw) 2713 2714 if len(missing_keys): 2715 warnings.warn( 2716 "Some non-standard WCS keywords were excluded: {} " 2717 "Use the ``relax`` kwarg to control this.".format( 2718 ', '.join(missing_keys)), 2719 AstropyWarning) 2720 # called when ``relax=None`` 2721 # This is different from the case of ``relax=False``. 2722 if any(self.wcs.ctype) and self.sip is not None: 2723 header = self._fix_ctype(header, add_sip=False, log_message=False) 2724 # Finally reset the key. This must be called after ``_fix_ctype``. 2725 if key is not None: 2726 self.wcs.alt = orig_key 2727 return header 2728 2729 def _fix_ctype(self, header, add_sip=True, log_message=True): 2730 """ 2731 Parameters 2732 ---------- 2733 header : `~astropy.io.fits.Header` 2734 FITS header. 2735 add_sip : bool 2736 Flag indicating whether "-SIP" should be added or removed from CTYPE keywords. 2737 2738 Remove "-SIP" from CTYPE when writing out a header with relax=False. 2739 This needs to be done outside ``to_header`` because ``to_header`` runs 2740 twice when ``relax=False`` and the second time ``relax`` is set to ``True`` 2741 to display the missing keywords. 2742 2743 If the user requested SIP distortion to be written out add "-SIP" to 2744 CTYPE if it is missing. 2745 """ 2746 2747 _add_sip_to_ctype = """ 2748 Inconsistent SIP distortion information is present in the current WCS: 2749 SIP coefficients were detected, but CTYPE is missing "-SIP" suffix, 2750 therefore the current WCS is internally inconsistent. 2751 2752 Because relax has been set to True, the resulting output WCS will have 2753 "-SIP" appended to CTYPE in order to make the header internally consistent. 2754 2755 However, this may produce incorrect astrometry in the output WCS, if 2756 in fact the current WCS is already distortion-corrected. 2757 2758 Therefore, if current WCS is already distortion-corrected (eg, drizzled) 2759 then SIP distortion components should not apply. In that case, for a WCS 2760 that is already distortion-corrected, please remove the SIP coefficients 2761 from the header. 2762 2763 """ 2764 if log_message: 2765 if add_sip: 2766 log.info(_add_sip_to_ctype) 2767 for i in range(1, self.naxis+1): 2768 # strip() must be called here to cover the case of alt key= " " 2769 kw = f'CTYPE{i}{self.wcs.alt}'.strip() 2770 if kw in header: 2771 if add_sip: 2772 val = header[kw].strip("-SIP") + "-SIP" 2773 else: 2774 val = header[kw].strip("-SIP") 2775 header[kw] = val 2776 else: 2777 continue 2778 return header 2779 2780 def to_header_string(self, relax=None): 2781 """ 2782 Identical to `to_header`, but returns a string containing the 2783 header cards. 2784 """ 2785 return str(self.to_header(relax)) 2786 2787 def footprint_to_file(self, filename='footprint.reg', color='green', 2788 width=2, coordsys=None): 2789 """ 2790 Writes out a `ds9`_ style regions file. It can be loaded 2791 directly by `ds9`_. 2792 2793 Parameters 2794 ---------- 2795 filename : str, optional 2796 Output file name - default is ``'footprint.reg'`` 2797 2798 color : str, optional 2799 Color to use when plotting the line. 2800 2801 width : int, optional 2802 Width of the region line. 2803 2804 coordsys : str, optional 2805 Coordinate system. If not specified (default), the ``radesys`` 2806 value is used. For all possible values, see 2807 http://ds9.si.edu/doc/ref/region.html#RegionFileFormat 2808 2809 """ 2810 comments = ('# Region file format: DS9 version 4.0 \n' 2811 '# global color=green font="helvetica 12 bold ' 2812 'select=1 highlite=1 edit=1 move=1 delete=1 ' 2813 'include=1 fixed=0 source\n') 2814 2815 coordsys = coordsys or self.wcs.radesys 2816 2817 if coordsys not in ('PHYSICAL', 'IMAGE', 'FK4', 'B1950', 'FK5', 2818 'J2000', 'GALACTIC', 'ECLIPTIC', 'ICRS', 'LINEAR', 2819 'AMPLIFIER', 'DETECTOR'): 2820 raise ValueError("Coordinate system '{}' is not supported. A valid" 2821 " one can be given with the 'coordsys' argument." 2822 .format(coordsys)) 2823 2824 with open(filename, mode='w') as f: 2825 f.write(comments) 2826 f.write(f'{coordsys}\n') 2827 f.write('polygon(') 2828 ftpr = self.calc_footprint() 2829 if ftpr is not None: 2830 ftpr.tofile(f, sep=',') 2831 f.write(f') # color={color}, width={width:d} \n') 2832 2833 def _get_naxis(self, header=None): 2834 _naxis = [] 2835 if (header is not None and 2836 not isinstance(header, (str, bytes))): 2837 for naxis in itertools.count(1): 2838 try: 2839 _naxis.append(header[f'NAXIS{naxis}']) 2840 except KeyError: 2841 break 2842 if len(_naxis) == 0: 2843 _naxis = [0, 0] 2844 elif len(_naxis) == 1: 2845 _naxis.append(0) 2846 self._naxis = _naxis 2847 2848 def printwcs(self): 2849 print(repr(self)) 2850 2851 def __repr__(self): 2852 ''' 2853 Return a short description. Simply porting the behavior from 2854 the `printwcs()` method. 2855 ''' 2856 description = ["WCS Keywords\n", 2857 f"Number of WCS axes: {self.naxis!r}"] 2858 sfmt = ' : ' + "".join(["{"+f"{i}"+"!r} " for i in range(self.naxis)]) 2859 2860 keywords = ['CTYPE', 'CRVAL', 'CRPIX'] 2861 values = [self.wcs.ctype, self.wcs.crval, self.wcs.crpix] 2862 for keyword, value in zip(keywords, values): 2863 description.append(keyword+sfmt.format(*value)) 2864 2865 if hasattr(self.wcs, 'pc'): 2866 for i in range(self.naxis): 2867 s = '' 2868 for j in range(self.naxis): 2869 s += ''.join(['PC', str(i+1), '_', str(j+1), ' ']) 2870 s += sfmt 2871 description.append(s.format(*self.wcs.pc[i])) 2872 s = 'CDELT' + sfmt 2873 description.append(s.format(*self.wcs.cdelt)) 2874 elif hasattr(self.wcs, 'cd'): 2875 for i in range(self.naxis): 2876 s = '' 2877 for j in range(self.naxis): 2878 s += "".join(['CD', str(i+1), '_', str(j+1), ' ']) 2879 s += sfmt 2880 description.append(s.format(*self.wcs.cd[i])) 2881 2882 description.append(f"NAXIS : {' '.join(map(str, self._naxis))}") 2883 return '\n'.join(description) 2884 2885 def get_axis_types(self): 2886 """ 2887 Similar to `self.wcsprm.axis_types <astropy.wcs.Wcsprm.axis_types>` 2888 but provides the information in a more Python-friendly format. 2889 2890 Returns 2891 ------- 2892 result : list of dict 2893 2894 Returns a list of dictionaries, one for each axis, each 2895 containing attributes about the type of that axis. 2896 2897 Each dictionary has the following keys: 2898 2899 - 'coordinate_type': 2900 2901 - None: Non-specific coordinate type. 2902 2903 - 'stokes': Stokes coordinate. 2904 2905 - 'celestial': Celestial coordinate (including ``CUBEFACE``). 2906 2907 - 'spectral': Spectral coordinate. 2908 2909 - 'scale': 2910 2911 - 'linear': Linear axis. 2912 2913 - 'quantized': Quantized axis (``STOKES``, ``CUBEFACE``). 2914 2915 - 'non-linear celestial': Non-linear celestial axis. 2916 2917 - 'non-linear spectral': Non-linear spectral axis. 2918 2919 - 'logarithmic': Logarithmic axis. 2920 2921 - 'tabular': Tabular axis. 2922 2923 - 'group' 2924 2925 - Group number, e.g. lookup table number 2926 2927 - 'number' 2928 2929 - For celestial axes: 2930 2931 - 0: Longitude coordinate. 2932 2933 - 1: Latitude coordinate. 2934 2935 - 2: ``CUBEFACE`` number. 2936 2937 - For lookup tables: 2938 2939 - the axis number in a multidimensional table. 2940 2941 ``CTYPEia`` in ``"4-3"`` form with unrecognized algorithm code will 2942 generate an error. 2943 """ 2944 if self.wcs is None: 2945 raise AttributeError( 2946 "This WCS object does not have a wcsprm object.") 2947 2948 coordinate_type_map = { 2949 0: None, 2950 1: 'stokes', 2951 2: 'celestial', 2952 3: 'spectral'} 2953 2954 scale_map = { 2955 0: 'linear', 2956 1: 'quantized', 2957 2: 'non-linear celestial', 2958 3: 'non-linear spectral', 2959 4: 'logarithmic', 2960 5: 'tabular'} 2961 2962 result = [] 2963 for axis_type in self.wcs.axis_types: 2964 subresult = {} 2965 2966 coordinate_type = (axis_type // 1000) % 10 2967 subresult['coordinate_type'] = coordinate_type_map[coordinate_type] 2968 2969 scale = (axis_type // 100) % 10 2970 subresult['scale'] = scale_map[scale] 2971 2972 group = (axis_type // 10) % 10 2973 subresult['group'] = group 2974 2975 number = axis_type % 10 2976 subresult['number'] = number 2977 2978 result.append(subresult) 2979 2980 return result 2981 2982 def __reduce__(self): 2983 """ 2984 Support pickling of WCS objects. This is done by serializing 2985 to an in-memory FITS file and dumping that as a string. 2986 """ 2987 2988 hdulist = self.to_fits(relax=True) 2989 2990 buffer = io.BytesIO() 2991 hdulist.writeto(buffer) 2992 2993 return (__WCS_unpickle__, 2994 (self.__class__, self.__dict__, buffer.getvalue(),)) 2995 2996 def dropaxis(self, dropax): 2997 """ 2998 Remove an axis from the WCS. 2999 3000 Parameters 3001 ---------- 3002 wcs : `~astropy.wcs.WCS` 3003 The WCS with naxis to be chopped to naxis-1 3004 dropax : int 3005 The index of the WCS to drop, counting from 0 (i.e., python convention, 3006 not FITS convention) 3007 3008 Returns 3009 ------- 3010 `~astropy.wcs.WCS` 3011 A new `~astropy.wcs.WCS` instance with one axis fewer 3012 """ 3013 inds = list(range(self.wcs.naxis)) 3014 inds.pop(dropax) 3015 3016 # axis 0 has special meaning to sub 3017 # if wcs.wcs.ctype == ['RA','DEC','VLSR'], you want 3018 # wcs.sub([1,2]) to get 'RA','DEC' back 3019 return self.sub([i+1 for i in inds]) 3020 3021 def swapaxes(self, ax0, ax1): 3022 """ 3023 Swap axes in a WCS. 3024 3025 Parameters 3026 ---------- 3027 wcs : `~astropy.wcs.WCS` 3028 The WCS to have its axes swapped 3029 ax0 : int 3030 ax1 : int 3031 The indices of the WCS to be swapped, counting from 0 (i.e., python 3032 convention, not FITS convention) 3033 3034 Returns 3035 ------- 3036 `~astropy.wcs.WCS` 3037 A new `~astropy.wcs.WCS` instance with the same number of axes, 3038 but two swapped 3039 """ 3040 inds = list(range(self.wcs.naxis)) 3041 inds[ax0], inds[ax1] = inds[ax1], inds[ax0] 3042 3043 return self.sub([i+1 for i in inds]) 3044 3045 def reorient_celestial_first(self): 3046 """ 3047 Reorient the WCS such that the celestial axes are first, followed by 3048 the spectral axis, followed by any others. 3049 Assumes at least celestial axes are present. 3050 """ 3051 return self.sub([WCSSUB_CELESTIAL, WCSSUB_SPECTRAL, WCSSUB_STOKES]) # Defined by C-ext # noqa: F821 E501 3052 3053 def slice(self, view, numpy_order=True): 3054 """ 3055 Slice a WCS instance using a Numpy slice. The order of the slice should 3056 be reversed (as for the data) compared to the natural WCS order. 3057 3058 Parameters 3059 ---------- 3060 view : tuple 3061 A tuple containing the same number of slices as the WCS system. 3062 The ``step`` method, the third argument to a slice, is not 3063 presently supported. 3064 numpy_order : bool 3065 Use numpy order, i.e. slice the WCS so that an identical slice 3066 applied to a numpy array will slice the array and WCS in the same 3067 way. If set to `False`, the WCS will be sliced in FITS order, 3068 meaning the first slice will be applied to the *last* numpy index 3069 but the *first* WCS axis. 3070 3071 Returns 3072 ------- 3073 wcs_new : `~astropy.wcs.WCS` 3074 A new resampled WCS axis 3075 """ 3076 if hasattr(view, '__len__') and len(view) > self.wcs.naxis: 3077 raise ValueError("Must have # of slices <= # of WCS axes") 3078 elif not hasattr(view, '__len__'): # view MUST be an iterable 3079 view = [view] 3080 3081 if not all(isinstance(x, slice) for x in view): 3082 # We need to drop some dimensions, but this may not always be 3083 # possible with .sub due to correlated axes, so instead we use the 3084 # generalized slicing infrastructure from astropy.wcs.wcsapi. 3085 return SlicedFITSWCS(self, view) 3086 3087 # NOTE: we could in principle use SlicedFITSWCS as above for all slicing, 3088 # but in the simple case where there are no axes dropped, we can just 3089 # create a full WCS object with updated WCS parameters which is faster 3090 # for this specific case and also backward-compatible. 3091 3092 wcs_new = self.deepcopy() 3093 if wcs_new.sip is not None: 3094 sip_crpix = wcs_new.sip.crpix.tolist() 3095 3096 for i, iview in enumerate(view): 3097 if iview.step is not None and iview.step < 0: 3098 raise NotImplementedError("Reversing an axis is not " 3099 "implemented.") 3100 3101 if numpy_order: 3102 wcs_index = self.wcs.naxis - 1 - i 3103 else: 3104 wcs_index = i 3105 3106 if iview.step is not None and iview.start is None: 3107 # Slice from "None" is equivalent to slice from 0 (but one 3108 # might want to downsample, so allow slices with 3109 # None,None,step or None,stop,step) 3110 iview = slice(0, iview.stop, iview.step) 3111 3112 if iview.start is not None: 3113 if iview.step not in (None, 1): 3114 crpix = self.wcs.crpix[wcs_index] 3115 cdelt = self.wcs.cdelt[wcs_index] 3116 # equivalently (keep this comment so you can compare eqns): 3117 # wcs_new.wcs.crpix[wcs_index] = 3118 # (crpix - iview.start)*iview.step + 0.5 - iview.step/2. 3119 crp = ((crpix - iview.start - 1.)/iview.step 3120 + 0.5 + 1./iview.step/2.) 3121 wcs_new.wcs.crpix[wcs_index] = crp 3122 if wcs_new.sip is not None: 3123 sip_crpix[wcs_index] = crp 3124 wcs_new.wcs.cdelt[wcs_index] = cdelt * iview.step 3125 else: 3126 wcs_new.wcs.crpix[wcs_index] -= iview.start 3127 if wcs_new.sip is not None: 3128 sip_crpix[wcs_index] -= iview.start 3129 3130 try: 3131 # range requires integers but the other attributes can also 3132 # handle arbitrary values, so this needs to be in a try/except. 3133 nitems = len(builtins.range(self._naxis[wcs_index])[iview]) 3134 except TypeError as exc: 3135 if 'indices must be integers' not in str(exc): 3136 raise 3137 warnings.warn("NAXIS{} attribute is not updated because at " 3138 "least one index ('{}') is no integer." 3139 "".format(wcs_index, iview), AstropyUserWarning) 3140 else: 3141 wcs_new._naxis[wcs_index] = nitems 3142 3143 if wcs_new.sip is not None: 3144 wcs_new.sip = Sip(self.sip.a, self.sip.b, self.sip.ap, self.sip.bp, 3145 sip_crpix) 3146 3147 return wcs_new 3148 3149 def __getitem__(self, item): 3150 # "getitem" is a shortcut for self.slice; it is very limited 3151 # there is no obvious and unambiguous interpretation of wcs[1,2,3] 3152 # We COULD allow wcs[1] to link to wcs.sub([2]) 3153 # (wcs[i] -> wcs.sub([i+1]) 3154 return self.slice(item) 3155 3156 def __iter__(self): 3157 # Having __getitem__ makes Python think WCS is iterable. However, 3158 # Python first checks whether __iter__ is present, so we can raise an 3159 # exception here. 3160 raise TypeError(f"'{self.__class__.__name__}' object is not iterable") 3161 3162 @property 3163 def axis_type_names(self): 3164 """ 3165 World names for each coordinate axis 3166 3167 Returns 3168 ------- 3169 list of str 3170 A list of names along each axis. 3171 """ 3172 names = list(self.wcs.cname) 3173 types = self.wcs.ctype 3174 for i in range(len(names)): 3175 if len(names[i]) > 0: 3176 continue 3177 names[i] = types[i].split('-')[0] 3178 return names 3179 3180 @property 3181 def celestial(self): 3182 """ 3183 A copy of the current WCS with only the celestial axes included 3184 """ 3185 return self.sub([WCSSUB_CELESTIAL]) # Defined by C-ext # noqa: F821 3186 3187 @property 3188 def is_celestial(self): 3189 return self.has_celestial and self.naxis == 2 3190 3191 @property 3192 def has_celestial(self): 3193 try: 3194 return self.wcs.lng >= 0 and self.wcs.lat >= 0 3195 except InconsistentAxisTypesError: 3196 return False 3197 3198 @property 3199 def spectral(self): 3200 """ 3201 A copy of the current WCS with only the spectral axes included 3202 """ 3203 return self.sub([WCSSUB_SPECTRAL]) # Defined by C-ext # noqa: F821 3204 3205 @property 3206 def is_spectral(self): 3207 return self.has_spectral and self.naxis == 1 3208 3209 @property 3210 def has_spectral(self): 3211 try: 3212 return self.wcs.spec >= 0 3213 except InconsistentAxisTypesError: 3214 return False 3215 3216 @property 3217 def has_distortion(self): 3218 """ 3219 Returns `True` if any distortion terms are present. 3220 """ 3221 return (self.sip is not None or 3222 self.cpdis1 is not None or self.cpdis2 is not None or 3223 self.det2im1 is not None and self.det2im2 is not None) 3224 3225 @property 3226 def pixel_scale_matrix(self): 3227 3228 try: 3229 cdelt = np.diag(self.wcs.get_cdelt()) 3230 pc = self.wcs.get_pc() 3231 except InconsistentAxisTypesError: 3232 try: 3233 # for non-celestial axes, get_cdelt doesn't work 3234 with warnings.catch_warnings(): 3235 warnings.filterwarnings( 3236 'ignore', 'cdelt will be ignored since cd is present', RuntimeWarning) 3237 cdelt = np.dot(self.wcs.cd, np.diag(self.wcs.cdelt)) 3238 except AttributeError: 3239 cdelt = np.diag(self.wcs.cdelt) 3240 3241 try: 3242 pc = self.wcs.pc 3243 except AttributeError: 3244 pc = 1 3245 3246 pccd = np.dot(cdelt, pc) 3247 3248 return pccd 3249 3250 def footprint_contains(self, coord, **kwargs): 3251 """ 3252 Determines if a given SkyCoord is contained in the wcs footprint. 3253 3254 Parameters 3255 ---------- 3256 coord : `~astropy.coordinates.SkyCoord` 3257 The coordinate to check if it is within the wcs coordinate. 3258 **kwargs : 3259 Additional arguments to pass to `~astropy.coordinates.SkyCoord.to_pixel` 3260 3261 Returns 3262 ------- 3263 response : bool 3264 True means the WCS footprint contains the coordinate, False means it does not. 3265 """ 3266 3267 return coord.contained_by(self, **kwargs) 3268 3269 3270def __WCS_unpickle__(cls, dct, fits_data): 3271 """ 3272 Unpickles a WCS object from a serialized FITS string. 3273 """ 3274 3275 self = cls.__new__(cls) 3276 self.__dict__.update(dct) 3277 3278 buffer = io.BytesIO(fits_data) 3279 hdulist = fits.open(buffer) 3280 3281 WCS.__init__(self, hdulist[0].header, hdulist) 3282 3283 return self 3284 3285 3286def find_all_wcs(header, relax=True, keysel=None, fix=True, 3287 translate_units='', 3288 _do_set=True): 3289 """ 3290 Find all the WCS transformations in the given header. 3291 3292 Parameters 3293 ---------- 3294 header : str or `~astropy.io.fits.Header` object. 3295 3296 relax : bool or int, optional 3297 Degree of permissiveness: 3298 3299 - `True` (default): Admit all recognized informal extensions of the 3300 WCS standard. 3301 3302 - `False`: Recognize only FITS keywords defined by the 3303 published WCS standard. 3304 3305 - `int`: a bit field selecting specific extensions to accept. 3306 See :ref:`astropy:relaxread` for details. 3307 3308 keysel : sequence of str, optional 3309 A list of flags used to select the keyword types considered by 3310 wcslib. When ``None``, only the standard image header 3311 keywords are considered (and the underlying wcspih() C 3312 function is called). To use binary table image array or pixel 3313 list keywords, *keysel* must be set. 3314 3315 Each element in the list should be one of the following strings: 3316 3317 - 'image': Image header keywords 3318 3319 - 'binary': Binary table image array keywords 3320 3321 - 'pixel': Pixel list keywords 3322 3323 Keywords such as ``EQUIna`` or ``RFRQna`` that are common to 3324 binary table image arrays and pixel lists (including 3325 ``WCSNna`` and ``TWCSna``) are selected by both 'binary' and 3326 'pixel'. 3327 3328 fix : bool, optional 3329 When `True` (default), call `~astropy.wcs.Wcsprm.fix` on 3330 the resulting objects to fix any non-standard uses in the 3331 header. `FITSFixedWarning` warnings will be emitted if any 3332 changes were made. 3333 3334 translate_units : str, optional 3335 Specify which potentially unsafe translations of non-standard 3336 unit strings to perform. By default, performs none. See 3337 `WCS.fix` for more information about this parameter. Only 3338 effective when ``fix`` is `True`. 3339 3340 Returns 3341 ------- 3342 wcses : list of `WCS` 3343 """ 3344 3345 if isinstance(header, (str, bytes)): 3346 header_string = header 3347 elif isinstance(header, fits.Header): 3348 header_string = header.tostring() 3349 else: 3350 raise TypeError( 3351 "header must be a string or astropy.io.fits.Header object") 3352 3353 keysel_flags = _parse_keysel(keysel) 3354 3355 if isinstance(header_string, str): 3356 header_bytes = header_string.encode('ascii') 3357 else: 3358 header_bytes = header_string 3359 3360 wcsprms = _wcs.find_all_wcs(header_bytes, relax, keysel_flags) 3361 3362 result = [] 3363 for wcsprm in wcsprms: 3364 subresult = WCS(fix=False, _do_set=False) 3365 subresult.wcs = wcsprm 3366 result.append(subresult) 3367 3368 if fix: 3369 subresult.fix(translate_units) 3370 3371 if _do_set: 3372 subresult.wcs.set() 3373 3374 return result 3375 3376 3377def validate(source): 3378 """ 3379 Prints a WCS validation report for the given FITS file. 3380 3381 Parameters 3382 ---------- 3383 source : str or file-like or `~astropy.io.fits.HDUList` 3384 The FITS file to validate. 3385 3386 Returns 3387 ------- 3388 results : list subclass instance 3389 The result is returned as nested lists. The first level 3390 corresponds to the HDUs in the given file. The next level has 3391 an entry for each WCS found in that header. The special 3392 subclass of list will pretty-print the results as a table when 3393 printed. 3394 3395 """ 3396 class _WcsValidateWcsResult(list): 3397 def __init__(self, key): 3398 self._key = key 3399 3400 def __repr__(self): 3401 result = [f" WCS key '{self._key or ' '}':"] 3402 if len(self): 3403 for entry in self: 3404 for i, line in enumerate(entry.splitlines()): 3405 if i == 0: 3406 initial_indent = ' - ' 3407 else: 3408 initial_indent = ' ' 3409 result.extend( 3410 textwrap.wrap( 3411 line, 3412 initial_indent=initial_indent, 3413 subsequent_indent=' ')) 3414 else: 3415 result.append(" No issues.") 3416 return '\n'.join(result) 3417 3418 class _WcsValidateHduResult(list): 3419 def __init__(self, hdu_index, hdu_name): 3420 self._hdu_index = hdu_index 3421 self._hdu_name = hdu_name 3422 list.__init__(self) 3423 3424 def __repr__(self): 3425 if len(self): 3426 if self._hdu_name: 3427 hdu_name = f' ({self._hdu_name})' 3428 else: 3429 hdu_name = '' 3430 result = [f'HDU {self._hdu_index}{hdu_name}:'] 3431 for wcs in self: 3432 result.append(repr(wcs)) 3433 return '\n'.join(result) 3434 return '' 3435 3436 class _WcsValidateResults(list): 3437 def __repr__(self): 3438 result = [] 3439 for hdu in self: 3440 content = repr(hdu) 3441 if len(content): 3442 result.append(content) 3443 return '\n\n'.join(result) 3444 3445 global __warningregistry__ 3446 3447 if isinstance(source, fits.HDUList): 3448 hdulist = source 3449 else: 3450 hdulist = fits.open(source) 3451 3452 results = _WcsValidateResults() 3453 3454 for i, hdu in enumerate(hdulist): 3455 hdu_results = _WcsValidateHduResult(i, hdu.name) 3456 results.append(hdu_results) 3457 3458 with warnings.catch_warnings(record=True) as warning_lines: 3459 wcses = find_all_wcs( 3460 hdu.header, relax=_wcs.WCSHDR_reject, 3461 fix=False, _do_set=False) 3462 3463 for wcs in wcses: 3464 wcs_results = _WcsValidateWcsResult(wcs.wcs.alt) 3465 hdu_results.append(wcs_results) 3466 3467 try: 3468 del __warningregistry__ 3469 except NameError: 3470 pass 3471 3472 with warnings.catch_warnings(record=True) as warning_lines: 3473 warnings.resetwarnings() 3474 warnings.simplefilter( 3475 "always", FITSFixedWarning, append=True) 3476 3477 try: 3478 WCS(hdu.header, 3479 key=wcs.wcs.alt or ' ', 3480 relax=_wcs.WCSHDR_reject, 3481 fix=True, _do_set=False) 3482 except WcsError as e: 3483 wcs_results.append(str(e)) 3484 3485 wcs_results.extend([str(x.message) for x in warning_lines]) 3486 3487 return results 3488