1# Licensed under a 3-clause BSD style license - see LICENSE.rst 2 3import io 4import os 5from datetime import datetime 6 7from packaging.version import Version 8import pytest 9import numpy as np 10from numpy.testing import ( 11 assert_allclose, assert_array_almost_equal, assert_array_almost_equal_nulp, 12 assert_array_equal) 13 14from astropy import wcs 15from astropy.wcs import _wcs # noqa 16from astropy import units as u 17from astropy.utils.data import ( 18 get_pkg_data_filenames, get_pkg_data_contents, get_pkg_data_filename) 19from astropy.utils.misc import NumpyRNGContext 20from astropy.utils.exceptions import ( 21 AstropyUserWarning, AstropyWarning, AstropyDeprecationWarning) 22from astropy.tests.helper import assert_quantity_allclose 23from astropy.io import fits 24from astropy.coordinates import SkyCoord 25from astropy.nddata import Cutout2D 26 27 28_WCSLIB_VER = Version(_wcs.__version__) 29 30 31# NOTE: User can choose to use system wcslib instead of bundled. 32def _check_v71_dateref_warnings(w, nmax=None): 33 if _WCSLIB_VER >= Version('7.1') and _WCSLIB_VER < Version('7.3') and w: 34 if nmax is None: 35 assert w 36 else: 37 assert len(w) == nmax 38 39 for item in w: 40 if (issubclass(item.category, wcs.FITSFixedWarning) and 41 str(item.message) == "'datfix' made the change " 42 "'Set DATE-REF to '1858-11-17' from MJD-REF'."): 43 break 44 else: 45 assert False, "No 'datfix' warning found" 46 47 48class TestMaps: 49 def setup(self): 50 # get the list of the hdr files that we want to test 51 self._file_list = list(get_pkg_data_filenames( 52 "data/maps", pattern="*.hdr")) 53 54 def test_consistency(self): 55 # Check to see that we actually have the list we expect, so that we 56 # do not get in a situation where the list is empty or incomplete and 57 # the tests still seem to pass correctly. 58 59 # how many do we expect to see? 60 n_data_files = 28 61 62 assert len(self._file_list) == n_data_files, ( 63 "test_spectra has wrong number data files: found {}, expected " 64 " {}".format(len(self._file_list), n_data_files)) 65 66 def test_maps(self): 67 for filename in self._file_list: 68 # use the base name of the file, so we get more useful messages 69 # for failing tests. 70 filename = os.path.basename(filename) 71 # Now find the associated file in the installed wcs test directory. 72 header = get_pkg_data_contents( 73 os.path.join("data", "maps", filename), encoding='binary') 74 # finally run the test. 75 wcsobj = wcs.WCS(header) 76 world = wcsobj.wcs_pix2world([[97, 97]], 1) 77 assert_array_almost_equal(world, [[285.0, -66.25]], decimal=1) 78 pix = wcsobj.wcs_world2pix([[285.0, -66.25]], 1) 79 assert_array_almost_equal(pix, [[97, 97]], decimal=0) 80 81 82class TestSpectra: 83 def setup(self): 84 self._file_list = list(get_pkg_data_filenames("data/spectra", 85 pattern="*.hdr")) 86 87 def test_consistency(self): 88 # Check to see that we actually have the list we expect, so that we 89 # do not get in a situation where the list is empty or incomplete and 90 # the tests still seem to pass correctly. 91 92 # how many do we expect to see? 93 n_data_files = 6 94 95 assert len(self._file_list) == n_data_files, ( 96 "test_spectra has wrong number data files: found {}, expected " 97 " {}".format(len(self._file_list), n_data_files)) 98 99 def test_spectra(self): 100 for filename in self._file_list: 101 # use the base name of the file, so we get more useful messages 102 # for failing tests. 103 filename = os.path.basename(filename) 104 # Now find the associated file in the installed wcs test directory. 105 header = get_pkg_data_contents( 106 os.path.join("data", "spectra", filename), encoding='binary') 107 # finally run the test. 108 with pytest.warns(None) as w: 109 all_wcs = wcs.find_all_wcs(header) 110 111 if _WCSLIB_VER >= Version('7.4'): 112 assert len(w) == 9 113 m = str(w.pop().message) 114 assert "'datfix' made the change 'Set MJD-OBS to 53925.853472 from DATE-OBS'." in m 115 else: 116 assert len(w) == 0 117 118 assert len(all_wcs) == 9 119 120 121def test_fixes(): 122 """ 123 From github issue #36 124 """ 125 header = get_pkg_data_contents('data/nonstandard_units.hdr', encoding='binary') 126 127 with pytest.raises(wcs.InvalidTransformError), pytest.warns(wcs.FITSFixedWarning) as w: 128 wcs.WCS(header, translate_units='dhs') 129 130 if Version('7.4') <=_WCSLIB_VER < Version('7.6'): 131 assert len(w) == 3 132 assert "'datfix' made the change 'Success'." in str(w.pop().message) 133 else: 134 assert len(w) == 2 135 136 first_wmsg = str(w[0].message) 137 assert 'unitfix' in first_wmsg and 'Hz' in first_wmsg and 'M/S' in first_wmsg 138 assert 'plane angle' in str(w[1].message) and 'm/s' in str(w[1].message) 139 140 141# Ignore "PV2_2 = 0.209028857410973 invalid keyvalue" warning seen on Windows. 142@pytest.mark.filterwarnings(r'ignore:PV2_2') 143def test_outside_sky(): 144 """ 145 From github issue #107 146 """ 147 header = get_pkg_data_contents( 148 'data/outside_sky.hdr', encoding='binary') 149 w = wcs.WCS(header) 150 151 assert np.all(np.isnan(w.wcs_pix2world([[100., 500.]], 0))) # outside sky 152 assert np.all(np.isnan(w.wcs_pix2world([[200., 200.]], 0))) # outside sky 153 assert not np.any(np.isnan(w.wcs_pix2world([[1000., 1000.]], 0))) 154 155 156def test_pix2world(): 157 """ 158 From github issue #1463 159 """ 160 # TODO: write this to test the expected output behavior of pix2world, 161 # currently this just makes sure it doesn't error out in unexpected ways 162 # (and compares `wcs.pc` and `result` values?) 163 filename = get_pkg_data_filename('data/sip2.fits') 164 with pytest.warns(wcs.FITSFixedWarning) as caught_warnings: 165 # this raises a warning unimportant for this testing the pix2world 166 # FITSFixedWarning(u'The WCS transformation has more axes (2) than 167 # the image it is associated with (0)') 168 ww = wcs.WCS(filename) 169 170 # might as well monitor for changing behavior 171 if Version('7.4') <=_WCSLIB_VER < Version('7.6'): 172 assert len(caught_warnings) == 2 173 else: 174 assert len(caught_warnings) == 1 175 176 n = 3 177 pixels = (np.arange(n) * np.ones((2, n))).T 178 result = ww.wcs_pix2world(pixels, 0, ra_dec_order=True) 179 180 # Catch #2791 181 ww.wcs_pix2world(pixels[..., 0], pixels[..., 1], 0, ra_dec_order=True) 182 183 # assuming that the data of sip2.fits doesn't change 184 answer = np.array([[0.00024976, 0.00023018], 185 [0.00023043, -0.00024997]]) 186 187 assert np.allclose(ww.wcs.pc, answer, atol=1.e-8) 188 189 answer = np.array([[202.39265216, 47.17756518], 190 [202.39335826, 47.17754619], 191 [202.39406436, 47.1775272]]) 192 193 assert np.allclose(result, answer, atol=1.e-8, rtol=1.e-10) 194 195 196def test_load_fits_path(): 197 fits_name = get_pkg_data_filename('data/sip.fits') 198 with pytest.warns(wcs.FITSFixedWarning): 199 wcs.WCS(fits_name) 200 201 202def test_dict_init(): 203 """ 204 Test that WCS can be initialized with a dict-like object 205 """ 206 207 # Dictionary with no actual WCS, returns identity transform 208 with pytest.warns(None) as wrng: 209 w = wcs.WCS({}) 210 _check_v71_dateref_warnings(wrng) 211 212 xp, yp = w.wcs_world2pix(41., 2., 1) 213 214 assert_array_almost_equal_nulp(xp, 41., 10) 215 assert_array_almost_equal_nulp(yp, 2., 10) 216 217 # Valid WCS 218 hdr = { 219 'CTYPE1': 'GLON-CAR', 220 'CTYPE2': 'GLAT-CAR', 221 'CUNIT1': 'deg', 222 'CUNIT2': 'deg', 223 'CRPIX1': 1, 224 'CRPIX2': 1, 225 'CRVAL1': 40., 226 'CRVAL2': 0., 227 'CDELT1': -0.1, 228 'CDELT2': 0.1 229 } 230 if _WCSLIB_VER >= Version('7.1'): 231 hdr['DATEREF'] = '1858-11-17' 232 233 with pytest.warns(None) as wrng: 234 w = wcs.WCS(hdr) 235 236 if _WCSLIB_VER >= Version('7.4'): 237 assert len(wrng) == 1 238 msg = str(wrng[0].message) 239 assert "'datfix' made the change 'Set MJDREF to 0.000000 from DATEREF'." in msg 240 else: 241 assert len(wrng) == 0 242 243 xp, yp = w.wcs_world2pix(41., 2., 0) 244 245 assert_array_almost_equal_nulp(xp, -10., 10) 246 assert_array_almost_equal_nulp(yp, 20., 10) 247 248 249def test_extra_kwarg(): 250 """ 251 Issue #444 252 """ 253 w = wcs.WCS() 254 with NumpyRNGContext(123456789): 255 data = np.random.rand(100, 2) 256 with pytest.raises(TypeError): 257 w.wcs_pix2world(data, origin=1) 258 259 260def test_3d_shapes(): 261 """ 262 Issue #444 263 """ 264 w = wcs.WCS(naxis=3) 265 with NumpyRNGContext(123456789): 266 data = np.random.rand(100, 3) 267 result = w.wcs_pix2world(data, 1) 268 assert result.shape == (100, 3) 269 result = w.wcs_pix2world( 270 data[..., 0], data[..., 1], data[..., 2], 1) 271 assert len(result) == 3 272 273 274def test_preserve_shape(): 275 w = wcs.WCS(naxis=2) 276 277 x = np.random.random((2, 3, 4)) 278 y = np.random.random((2, 3, 4)) 279 280 xw, yw = w.wcs_pix2world(x, y, 1) 281 282 assert xw.shape == (2, 3, 4) 283 assert yw.shape == (2, 3, 4) 284 285 xp, yp = w.wcs_world2pix(x, y, 1) 286 287 assert xp.shape == (2, 3, 4) 288 assert yp.shape == (2, 3, 4) 289 290 291def test_broadcasting(): 292 w = wcs.WCS(naxis=2) 293 294 x = np.random.random((2, 3, 4)) 295 y = 1 296 297 xp, yp = w.wcs_world2pix(x, y, 1) 298 299 assert xp.shape == (2, 3, 4) 300 assert yp.shape == (2, 3, 4) 301 302 303def test_shape_mismatch(): 304 w = wcs.WCS(naxis=2) 305 306 x = np.random.random((2, 3, 4)) 307 y = np.random.random((3, 2, 4)) 308 309 with pytest.raises(ValueError) as exc: 310 xw, yw = w.wcs_pix2world(x, y, 1) 311 assert exc.value.args[0] == "Coordinate arrays are not broadcastable to each other" 312 313 with pytest.raises(ValueError) as exc: 314 xp, yp = w.wcs_world2pix(x, y, 1) 315 assert exc.value.args[0] == "Coordinate arrays are not broadcastable to each other" 316 317 # There are some ambiguities that need to be worked around when 318 # naxis == 1 319 w = wcs.WCS(naxis=1) 320 321 x = np.random.random((42, 1)) 322 xw = w.wcs_pix2world(x, 1) 323 assert xw.shape == (42, 1) 324 325 x = np.random.random((42,)) 326 xw, = w.wcs_pix2world(x, 1) 327 assert xw.shape == (42,) 328 329 330def test_invalid_shape(): 331 # Issue #1395 332 w = wcs.WCS(naxis=2) 333 334 xy = np.random.random((2, 3)) 335 with pytest.raises(ValueError) as exc: 336 w.wcs_pix2world(xy, 1) 337 assert exc.value.args[0] == 'When providing two arguments, the array must be of shape (N, 2)' 338 339 xy = np.random.random((2, 1)) 340 with pytest.raises(ValueError) as exc: 341 w.wcs_pix2world(xy, 1) 342 assert exc.value.args[0] == 'When providing two arguments, the array must be of shape (N, 2)' 343 344 345def test_warning_about_defunct_keywords(): 346 header = get_pkg_data_contents('data/defunct_keywords.hdr', encoding='binary') 347 if Version('7.4') <=_WCSLIB_VER < Version('7.6'): 348 n_warn = 5 349 else: 350 n_warn = 4 351 352 # Make sure the warnings come out every time... 353 for _ in range(2): 354 with pytest.warns(wcs.FITSFixedWarning) as w: 355 wcs.WCS(header) 356 357 assert len(w) == n_warn 358 # 7.4 adds a fifth warning "'datfix' made the change 'Success'." 359 for item in w[:4]: 360 assert 'PCi_ja' in str(item.message) 361 362 363def test_warning_about_defunct_keywords_exception(): 364 header = get_pkg_data_contents('data/defunct_keywords.hdr', encoding='binary') 365 with pytest.warns(wcs.FITSFixedWarning): 366 wcs.WCS(header) 367 368 369def test_to_header_string(): 370 hdrstr = ( 371 "WCSAXES = 2 / Number of coordinate axes ", 372 "CRPIX1 = 0.0 / Pixel coordinate of reference point ", 373 "CRPIX2 = 0.0 / Pixel coordinate of reference point ", 374 "CDELT1 = 1.0 / Coordinate increment at reference point ", 375 "CDELT2 = 1.0 / Coordinate increment at reference point ", 376 "CRVAL1 = 0.0 / Coordinate value at reference point ", 377 "CRVAL2 = 0.0 / Coordinate value at reference point ", 378 "LATPOLE = 90.0 / [deg] Native latitude of celestial pole ", 379 ) 380 381 if _WCSLIB_VER >= Version('7.3'): 382 hdrstr += ( 383 "MJDREF = 0.0 / [d] MJD of fiducial time ", 384 ) 385 386 elif _WCSLIB_VER >= Version('7.1'): 387 hdrstr += ( 388 "DATEREF = '1858-11-17' / ISO-8601 fiducial time ", 389 "MJDREFI = 0.0 / [d] MJD of fiducial time, integer part ", 390 "MJDREFF = 0.0 / [d] MJD of fiducial time, fractional part " 391 ) 392 393 hdrstr += ("END", ) 394 395 header_string = ''.join(hdrstr) 396 397 w = wcs.WCS() 398 h0 = fits.Header.fromstring(w.to_header_string().strip()) 399 if 'COMMENT' in h0: 400 del h0['COMMENT'] 401 if '' in h0: 402 del h0[''] 403 h1 = fits.Header.fromstring(header_string.strip()) 404 assert dict(h0) == dict(h1) 405 406 407def test_to_fits(): 408 nrec = 11 if _WCSLIB_VER >= Version('7.1') else 8 409 if _WCSLIB_VER < Version('7.1'): 410 nrec = 8 411 elif _WCSLIB_VER < Version('7.3'): 412 nrec = 11 413 else: 414 nrec = 9 415 416 w = wcs.WCS() 417 header_string = w.to_header() 418 wfits = w.to_fits() 419 assert isinstance(wfits, fits.HDUList) 420 assert isinstance(wfits[0], fits.PrimaryHDU) 421 assert header_string == wfits[0].header[-nrec:] 422 423 424def test_to_header_warning(): 425 fits_name = get_pkg_data_filename('data/sip.fits') 426 with pytest.warns(wcs.FITSFixedWarning): 427 x = wcs.WCS(fits_name) 428 with pytest.warns(AstropyWarning, match='A_ORDER') as w: 429 x.to_header() 430 assert len(w) == 1 431 432 433def test_no_comments_in_header(): 434 w = wcs.WCS() 435 header = w.to_header() 436 assert w.wcs.alt not in header 437 assert 'COMMENT' + w.wcs.alt.strip() not in header 438 assert 'COMMENT' not in header 439 wkey = 'P' 440 header = w.to_header(key=wkey) 441 assert wkey not in header 442 assert 'COMMENT' not in header 443 assert 'COMMENT' + w.wcs.alt.strip() not in header 444 445 446def test_find_all_wcs_crash(): 447 """ 448 Causes a double free without a recent fix in wcslib_wrap.C 449 """ 450 with open(get_pkg_data_filename("data/too_many_pv.hdr")) as fd: 451 header = fd.read() 452 # We have to set fix=False here, because one of the fixing tasks is to 453 # remove redundant SCAMP distortion parameters when SIP distortion 454 # parameters are also present. 455 with pytest.raises(wcs.InvalidTransformError), pytest.warns(wcs.FITSFixedWarning): 456 wcs.find_all_wcs(header, fix=False) 457 458 459# NOTE: Warning bubbles up from C layer during wcs.validate() and 460# is hard to catch, so we just ignore it. 461@pytest.mark.filterwarnings("ignore") 462def test_validate(): 463 results = wcs.validate(get_pkg_data_filename("data/validate.fits")) 464 results_txt = sorted(set([x.strip() for x in repr(results).splitlines()])) 465 if _WCSLIB_VER >= Version('7.6'): 466 filename = 'data/validate.7.6.txt' 467 elif _WCSLIB_VER >= Version('7.4'): 468 filename = 'data/validate.7.4.txt' 469 elif _WCSLIB_VER >= Version('6.0'): 470 filename = 'data/validate.6.txt' 471 elif _WCSLIB_VER >= Version('5.13'): 472 filename = 'data/validate.5.13.txt' 473 elif _WCSLIB_VER >= Version('5.0'): 474 filename = 'data/validate.5.0.txt' 475 else: 476 filename = 'data/validate.txt' 477 with open(get_pkg_data_filename(filename), "r") as fd: 478 lines = fd.readlines() 479 assert sorted(set([x.strip() for x in lines])) == results_txt 480 481 482def test_validate_with_2_wcses(): 483 # From Issue #2053 484 with pytest.warns(AstropyUserWarning): 485 results = wcs.validate(get_pkg_data_filename("data/2wcses.hdr")) 486 487 assert "WCS key 'A':" in str(results) 488 489 490def test_crpix_maps_to_crval(): 491 twcs = wcs.WCS(naxis=2) 492 twcs.wcs.crval = [251.29, 57.58] 493 twcs.wcs.cdelt = [1, 1] 494 twcs.wcs.crpix = [507, 507] 495 twcs.wcs.pc = np.array([[7.7e-6, 3.3e-5], [3.7e-5, -6.8e-6]]) 496 twcs._naxis = [1014, 1014] 497 twcs.wcs.ctype = ['RA---TAN-SIP', 'DEC--TAN-SIP'] 498 a = np.array( 499 [[0, 0, 5.33092692e-08, 3.73753773e-11, -2.02111473e-13], 500 [0, 2.44084308e-05, 2.81394789e-11, 5.17856895e-13, 0.0], 501 [-2.41334657e-07, 1.29289255e-10, 2.35753629e-14, 0.0, 0.0], 502 [-2.37162007e-10, 5.43714947e-13, 0.0, 0.0, 0.0], 503 [-2.81029767e-13, 0.0, 0.0, 0.0, 0.0]] 504 ) 505 b = np.array( 506 [[0, 0, 2.99270374e-05, -2.38136074e-10, 7.23205168e-13], 507 [0, -1.71073858e-07, 6.31243431e-11, -5.16744347e-14, 0.0], 508 [6.95458963e-06, -3.08278961e-10, -1.75800917e-13, 0.0, 0.0], 509 [3.51974159e-11, 5.60993016e-14, 0.0, 0.0, 0.0], 510 [-5.92438525e-13, 0.0, 0.0, 0.0, 0.0]] 511 ) 512 twcs.sip = wcs.Sip(a, b, None, None, twcs.wcs.crpix) 513 twcs.wcs.set() 514 pscale = np.sqrt(wcs.utils.proj_plane_pixel_area(twcs)) 515 516 # test that CRPIX maps to CRVAL: 517 assert_allclose( 518 twcs.wcs_pix2world(*twcs.wcs.crpix, 1), twcs.wcs.crval, 519 rtol=0.0, atol=1e-6 * pscale 520 ) 521 522 # test that CRPIX maps to CRVAL: 523 assert_allclose( 524 twcs.all_pix2world(*twcs.wcs.crpix, 1), twcs.wcs.crval, 525 rtol=0.0, atol=1e-6 * pscale 526 ) 527 528 529def test_all_world2pix(fname=None, ext=0, 530 tolerance=1.0e-4, origin=0, 531 random_npts=25000, 532 adaptive=False, maxiter=20, 533 detect_divergence=True): 534 """Test all_world2pix, iterative inverse of all_pix2world""" 535 536 # Open test FITS file: 537 if fname is None: 538 fname = get_pkg_data_filename('data/j94f05bgq_flt.fits') 539 ext = ('SCI', 1) 540 if not os.path.isfile(fname): 541 raise OSError(f"Input file '{fname:s}' to 'test_all_world2pix' not found.") 542 h = fits.open(fname) 543 w = wcs.WCS(h[ext].header, h) 544 h.close() 545 del h 546 547 crpix = w.wcs.crpix 548 ncoord = crpix.shape[0] 549 550 # Assume that CRPIX is at the center of the image and that the image has 551 # a power-of-2 number of pixels along each axis. Only use the central 552 # 1/64 for this testing purpose: 553 naxesi_l = list((7. / 16 * crpix).astype(int)) 554 naxesi_u = list((9. / 16 * crpix).astype(int)) 555 556 # Generate integer indices of pixels (image grid): 557 img_pix = np.dstack([i.flatten() for i in 558 np.meshgrid(*map(range, naxesi_l, naxesi_u))])[0] 559 560 # Generage random data (in image coordinates): 561 with NumpyRNGContext(123456789): 562 rnd_pix = np.random.rand(random_npts, ncoord) 563 564 # Scale random data to cover the central part of the image 565 mwidth = 2 * (crpix * 1. / 8) 566 rnd_pix = crpix - 0.5 * mwidth + (mwidth - 1) * rnd_pix 567 568 # Reference pixel coordinates in image coordinate system (CS): 569 test_pix = np.append(img_pix, rnd_pix, axis=0) 570 # Reference pixel coordinates in sky CS using forward transformation: 571 all_world = w.all_pix2world(test_pix, origin) 572 573 try: 574 runtime_begin = datetime.now() 575 # Apply the inverse iterative process to pixels in world coordinates 576 # to recover the pixel coordinates in image space. 577 all_pix = w.all_world2pix( 578 all_world, origin, tolerance=tolerance, adaptive=adaptive, 579 maxiter=maxiter, detect_divergence=detect_divergence) 580 runtime_end = datetime.now() 581 except wcs.wcs.NoConvergence as e: 582 runtime_end = datetime.now() 583 ndiv = 0 584 if e.divergent is not None: 585 ndiv = e.divergent.shape[0] 586 print(f"There are {ndiv} diverging solutions.") 587 print(f"Indices of diverging solutions:\n{e.divergent}") 588 print(f"Diverging solutions:\n{e.best_solution[e.divergent]}\n") 589 print("Mean radius of the diverging solutions: {}" 590 .format(np.mean( 591 np.linalg.norm(e.best_solution[e.divergent], axis=1)))) 592 print("Mean accuracy of the diverging solutions: {}\n" 593 .format(np.mean( 594 np.linalg.norm(e.accuracy[e.divergent], axis=1)))) 595 else: 596 print("There are no diverging solutions.") 597 598 nslow = 0 599 if e.slow_conv is not None: 600 nslow = e.slow_conv.shape[0] 601 print(f"There are {nslow} slowly converging solutions.") 602 print(f"Indices of slowly converging solutions:\n{e.slow_conv}") 603 print(f"Slowly converging solutions:\n{e.best_solution[e.slow_conv]}\n") 604 else: 605 print("There are no slowly converging solutions.\n") 606 607 print("There are {} converged solutions." 608 .format(e.best_solution.shape[0] - ndiv - nslow)) 609 print(f"Best solutions (all points):\n{e.best_solution}") 610 print(f"Accuracy:\n{e.accuracy}\n") 611 print("\nFinished running 'test_all_world2pix' with errors.\n" 612 "ERROR: {}\nRun time: {}\n" 613 .format(e.args[0], runtime_end - runtime_begin)) 614 raise e 615 616 # Compute differences between reference pixel coordinates and 617 # pixel coordinates (in image space) recovered from reference 618 # pixels in world coordinates: 619 errors = np.sqrt(np.sum(np.power(all_pix - test_pix, 2), axis=1)) 620 meanerr = np.mean(errors) 621 maxerr = np.amax(errors) 622 print("\nFinished running 'test_all_world2pix'.\n" 623 "Mean error = {:e} (Max error = {:e})\n" 624 "Run time: {}\n" 625 .format(meanerr, maxerr, runtime_end - runtime_begin)) 626 627 assert(maxerr < 2.0 * tolerance) 628 629 630def test_scamp_sip_distortion_parameters(): 631 """ 632 Test parsing of WCS parameters with redundant SIP and SCAMP distortion 633 parameters. 634 """ 635 header = get_pkg_data_contents('data/validate.fits', encoding='binary') 636 with pytest.warns(wcs.FITSFixedWarning): 637 w = wcs.WCS(header) 638 # Just check that this doesn't raise an exception. 639 w.all_pix2world(0, 0, 0) 640 641 642def test_fixes2(): 643 """ 644 From github issue #1854 645 """ 646 header = get_pkg_data_contents( 647 'data/nonstandard_units.hdr', encoding='binary') 648 with pytest.raises(wcs.InvalidTransformError): 649 wcs.WCS(header, fix=False) 650 651 652def test_unit_normalization(): 653 """ 654 From github issue #1918 655 """ 656 header = get_pkg_data_contents( 657 'data/unit.hdr', encoding='binary') 658 w = wcs.WCS(header) 659 assert w.wcs.cunit[2] == 'm/s' 660 661 662def test_footprint_to_file(tmpdir): 663 """ 664 From github issue #1912 665 """ 666 # Arbitrary keywords from real data 667 hdr = {'CTYPE1': 'RA---ZPN', 'CRUNIT1': 'deg', 668 'CRPIX1': -3.3495999e+02, 'CRVAL1': 3.185790700000e+02, 669 'CTYPE2': 'DEC--ZPN', 'CRUNIT2': 'deg', 670 'CRPIX2': 3.0453999e+03, 'CRVAL2': 4.388538000000e+01, 671 'PV2_1': 1., 'PV2_3': 220., 'NAXIS1': 2048, 'NAXIS2': 1024} 672 w = wcs.WCS(hdr) 673 674 testfile = str(tmpdir.join('test.txt')) 675 w.footprint_to_file(testfile) 676 677 with open(testfile, 'r') as f: 678 lines = f.readlines() 679 680 assert len(lines) == 4 681 assert lines[2] == 'ICRS\n' 682 assert 'color=green' in lines[3] 683 684 w.footprint_to_file(testfile, coordsys='FK5', color='red') 685 686 with open(testfile, 'r') as f: 687 lines = f.readlines() 688 689 assert len(lines) == 4 690 assert lines[2] == 'FK5\n' 691 assert 'color=red' in lines[3] 692 693 with pytest.raises(ValueError): 694 w.footprint_to_file(testfile, coordsys='FOO') 695 696 del hdr['NAXIS1'] 697 del hdr['NAXIS2'] 698 w = wcs.WCS(hdr) 699 with pytest.warns(AstropyUserWarning): 700 w.footprint_to_file(testfile) 701 702 703# Ignore FITSFixedWarning about keyrecords following the END keyrecord were 704# ignored, which comes from src/astropy_wcs.c . Only a blind catch like this 705# seems to work when pytest warnings are turned into exceptions. 706@pytest.mark.filterwarnings('ignore') 707def test_validate_faulty_wcs(): 708 """ 709 From github issue #2053 710 """ 711 h = fits.Header() 712 # Illegal WCS: 713 h['RADESYSA'] = 'ICRS' 714 h['PV2_1'] = 1.0 715 hdu = fits.PrimaryHDU([[0]], header=h) 716 hdulist = fits.HDUList([hdu]) 717 # Check that this doesn't raise a NameError exception 718 wcs.validate(hdulist) 719 720 721def test_error_message(): 722 header = get_pkg_data_contents( 723 'data/invalid_header.hdr', encoding='binary') 724 725 with pytest.raises(wcs.InvalidTransformError): 726 # Both lines are in here, because 0.4 calls .set within WCS.__init__, 727 # whereas 0.3 and earlier did not. 728 with pytest.warns(wcs.FITSFixedWarning): 729 w = wcs.WCS(header, _do_set=False) 730 w.all_pix2world([[536.0, 894.0]], 0) 731 732 733def test_out_of_bounds(): 734 # See #2107 735 header = get_pkg_data_contents('data/zpn-hole.hdr', encoding='binary') 736 w = wcs.WCS(header) 737 738 ra, dec = w.wcs_pix2world(110, 110, 0) 739 740 assert np.isnan(ra) 741 assert np.isnan(dec) 742 743 ra, dec = w.wcs_pix2world(0, 0, 0) 744 745 assert not np.isnan(ra) 746 assert not np.isnan(dec) 747 748 749def test_calc_footprint_1(): 750 fits = get_pkg_data_filename('data/sip.fits') 751 with pytest.warns(wcs.FITSFixedWarning): 752 w = wcs.WCS(fits) 753 754 axes = (1000, 1051) 755 ref = np.array([[202.39314493, 47.17753352], 756 [202.71885939, 46.94630488], 757 [202.94631893, 47.15855022], 758 [202.72053428, 47.37893142]]) 759 footprint = w.calc_footprint(axes=axes) 760 assert_allclose(footprint, ref) 761 762 763def test_calc_footprint_2(): 764 """ Test calc_footprint without distortion. """ 765 fits = get_pkg_data_filename('data/sip.fits') 766 with pytest.warns(wcs.FITSFixedWarning): 767 w = wcs.WCS(fits) 768 769 axes = (1000, 1051) 770 ref = np.array([[202.39265216, 47.17756518], 771 [202.7469062, 46.91483312], 772 [203.11487481, 47.14359319], 773 [202.76092671, 47.40745948]]) 774 footprint = w.calc_footprint(axes=axes, undistort=False) 775 assert_allclose(footprint, ref) 776 777 778def test_calc_footprint_3(): 779 """ Test calc_footprint with corner of the pixel.""" 780 w = wcs.WCS() 781 w.wcs.ctype = ["GLON-CAR", "GLAT-CAR"] 782 w.wcs.crpix = [1.5, 5.5] 783 w.wcs.cdelt = [-0.1, 0.1] 784 axes = (2, 10) 785 ref = np.array([[0.1, -0.5], 786 [0.1, 0.5], 787 [359.9, 0.5], 788 [359.9, -0.5]]) 789 790 footprint = w.calc_footprint(axes=axes, undistort=False, center=False) 791 assert_allclose(footprint, ref) 792 793 794def test_sip(): 795 # See #2107 796 header = get_pkg_data_contents('data/irac_sip.hdr', encoding='binary') 797 w = wcs.WCS(header) 798 799 x0, y0 = w.sip_pix2foc(200, 200, 0) 800 801 assert_allclose(72, x0, 1e-3) 802 assert_allclose(72, y0, 1e-3) 803 804 x1, y1 = w.sip_foc2pix(x0, y0, 0) 805 806 assert_allclose(200, x1, 1e-3) 807 assert_allclose(200, y1, 1e-3) 808 809 810def test_sub_3d_with_sip(): 811 # See #10527 812 header = get_pkg_data_contents('data/irac_sip.hdr', encoding='binary') 813 header = fits.Header.fromstring(header) 814 header['NAXIS'] = 3 815 header.set('NAXIS3', 64, after=header.index('NAXIS2')) 816 w = wcs.WCS(header, naxis=2) 817 assert w.naxis == 2 818 819 820def test_printwcs(capsys): 821 """ 822 Just make sure that it runs 823 """ 824 h = get_pkg_data_contents( 825 'data/spectra/orion-freq-1.hdr', encoding='binary') 826 with pytest.warns(wcs.FITSFixedWarning): 827 w = wcs.WCS(h) 828 w.printwcs() 829 captured = capsys.readouterr() 830 assert 'WCS Keywords' in captured.out 831 h = get_pkg_data_contents('data/3d_cd.hdr', encoding='binary') 832 w = wcs.WCS(h) 833 w.printwcs() 834 captured = capsys.readouterr() 835 assert 'WCS Keywords' in captured.out 836 837 838def test_invalid_spherical(): 839 header = """ 840SIMPLE = T / conforms to FITS standard 841BITPIX = 8 / array data type 842WCSAXES = 2 / no comment 843CTYPE1 = 'RA---TAN' / TAN (gnomic) projection 844CTYPE2 = 'DEC--TAN' / TAN (gnomic) projection 845EQUINOX = 2000.0 / Equatorial coordinates definition (yr) 846LONPOLE = 180.0 / no comment 847LATPOLE = 0.0 / no comment 848CRVAL1 = 16.0531567459 / RA of reference point 849CRVAL2 = 23.1148929108 / DEC of reference point 850CRPIX1 = 2129 / X reference pixel 851CRPIX2 = 1417 / Y reference pixel 852CUNIT1 = 'deg ' / X pixel scale units 853CUNIT2 = 'deg ' / Y pixel scale units 854CD1_1 = -0.00912247310646 / Transformation matrix 855CD1_2 = -0.00250608809647 / no comment 856CD2_1 = 0.00250608809647 / no comment 857CD2_2 = -0.00912247310646 / no comment 858IMAGEW = 4256 / Image width, in pixels. 859IMAGEH = 2832 / Image height, in pixels. 860 """ 861 862 f = io.StringIO(header) 863 header = fits.Header.fromtextfile(f) 864 865 w = wcs.WCS(header) 866 x, y = w.wcs_world2pix(211, -26, 0) 867 assert np.isnan(x) and np.isnan(y) 868 869 870def test_no_iteration(): 871 872 # Regression test for #3066 873 874 w = wcs.WCS(naxis=2) 875 876 with pytest.raises(TypeError) as exc: 877 iter(w) 878 assert exc.value.args[0] == "'WCS' object is not iterable" 879 880 class NewWCS(wcs.WCS): 881 pass 882 883 w = NewWCS(naxis=2) 884 885 with pytest.raises(TypeError) as exc: 886 iter(w) 887 assert exc.value.args[0] == "'NewWCS' object is not iterable" 888 889 890@pytest.mark.skipif('_wcs.__version__[0] < "5"', 891 reason="TPV only works with wcslib 5.x or later") 892def test_sip_tpv_agreement(): 893 sip_header = get_pkg_data_contents( 894 os.path.join("data", "siponly.hdr"), encoding='binary') 895 tpv_header = get_pkg_data_contents( 896 os.path.join("data", "tpvonly.hdr"), encoding='binary') 897 898 with pytest.warns(wcs.FITSFixedWarning): 899 w_sip = wcs.WCS(sip_header) 900 w_tpv = wcs.WCS(tpv_header) 901 902 assert_array_almost_equal( 903 w_sip.all_pix2world([w_sip.wcs.crpix], 1), 904 w_tpv.all_pix2world([w_tpv.wcs.crpix], 1)) 905 906 w_sip2 = wcs.WCS(w_sip.to_header()) 907 w_tpv2 = wcs.WCS(w_tpv.to_header()) 908 909 assert_array_almost_equal( 910 w_sip.all_pix2world([w_sip.wcs.crpix], 1), 911 w_sip2.all_pix2world([w_sip.wcs.crpix], 1)) 912 assert_array_almost_equal( 913 w_tpv.all_pix2world([w_sip.wcs.crpix], 1), 914 w_tpv2.all_pix2world([w_sip.wcs.crpix], 1)) 915 assert_array_almost_equal( 916 w_sip2.all_pix2world([w_sip.wcs.crpix], 1), 917 w_tpv2.all_pix2world([w_tpv.wcs.crpix], 1)) 918 919 920@pytest.mark.skipif('_wcs.__version__[0] < "5"', 921 reason="TPV only works with wcslib 5.x or later") 922def test_tpv_copy(): 923 # See #3904 924 925 tpv_header = get_pkg_data_contents( 926 os.path.join("data", "tpvonly.hdr"), encoding='binary') 927 928 with pytest.warns(wcs.FITSFixedWarning): 929 w_tpv = wcs.WCS(tpv_header) 930 931 ra, dec = w_tpv.wcs_pix2world([0, 100, 200], [0, -100, 200], 0) 932 assert ra[0] != ra[1] and ra[1] != ra[2] 933 assert dec[0] != dec[1] and dec[1] != dec[2] 934 935 936def test_hst_wcs(): 937 path = get_pkg_data_filename("data/dist_lookup.fits.gz") 938 939 with fits.open(path) as hdulist: 940 # wcslib will complain about the distortion parameters if they 941 # weren't correctly deleted from the header 942 w = wcs.WCS(hdulist[1].header, hdulist) 943 944 # Check pixel scale and area 945 assert_quantity_allclose( 946 w.proj_plane_pixel_scales(), [1.38484378e-05, 1.39758488e-05] * u.deg) 947 assert_quantity_allclose( 948 w.proj_plane_pixel_area(), 1.93085492e-10 * (u.deg * u.deg)) 949 950 # Exercise the main transformation functions, mainly just for 951 # coverage 952 w.p4_pix2foc([0, 100, 200], [0, -100, 200], 0) 953 w.det2im([0, 100, 200], [0, -100, 200], 0) 954 955 w.cpdis1 = w.cpdis1 956 w.cpdis2 = w.cpdis2 957 958 w.det2im1 = w.det2im1 959 w.det2im2 = w.det2im2 960 961 w.sip = w.sip 962 963 w.cpdis1.cdelt = w.cpdis1.cdelt 964 w.cpdis1.crpix = w.cpdis1.crpix 965 w.cpdis1.crval = w.cpdis1.crval 966 w.cpdis1.data = w.cpdis1.data 967 968 assert w.sip.a_order == 4 969 assert w.sip.b_order == 4 970 assert w.sip.ap_order == 0 971 assert w.sip.bp_order == 0 972 assert_array_equal(w.sip.crpix, [2048., 1024.]) 973 wcs.WCS(hdulist[1].header, hdulist) 974 975 976def test_cpdis_comments(): 977 path = get_pkg_data_filename("data/dist_lookup.fits.gz") 978 979 f = fits.open(path) 980 w = wcs.WCS(f[1].header, f) 981 hdr = w.to_fits()[0].header 982 f.close() 983 984 wcscards = list(hdr['CPDIS*'].cards) + list(hdr['DP*'].cards) 985 wcsdict = {k: (v, c) for k, v, c in wcscards} 986 987 refcards = [ 988 ('CPDIS1', 'LOOKUP', 'Prior distortion function type'), 989 ('DP1.EXTVER', 1.0, 'Version number of WCSDVARR extension'), 990 ('DP1.NAXES', 2.0, 'Number of independent variables in CPDIS function'), 991 ('DP1.AXIS.1', 1.0, 'Axis number of the 1st variable in a CPDIS function'), 992 ('DP1.AXIS.2', 2.0, 'Axis number of the 2nd variable in a CPDIS function'), 993 ('CPDIS2', 'LOOKUP', 'Prior distortion function type'), 994 ('DP2.EXTVER', 2.0, 'Version number of WCSDVARR extension'), 995 ('DP2.NAXES', 2.0, 'Number of independent variables in CPDIS function'), 996 ('DP2.AXIS.1', 1.0, 'Axis number of the 1st variable in a CPDIS function'), 997 ('DP2.AXIS.2', 2.0, 'Axis number of the 2nd variable in a CPDIS function'), 998 ] 999 1000 assert len(wcsdict) == len(refcards) 1001 1002 for k, v, c in refcards: 1003 assert wcsdict[k] == (v, c) 1004 1005 1006def test_d2im_comments(): 1007 path = get_pkg_data_filename("data/ie6d07ujq_wcs.fits") 1008 1009 f = fits.open(path) 1010 with pytest.warns(wcs.FITSFixedWarning): 1011 w = wcs.WCS(f[0].header, f) 1012 f.close() 1013 wcscards = list(w.to_fits()[0].header['D2IM*'].cards) 1014 wcsdict = {k: (v, c) for k, v, c in wcscards} 1015 1016 refcards = [ 1017 ('D2IMDIS1', 'LOOKUP', 'Detector to image correction type'), 1018 ('D2IM1.EXTVER', 1.0, 'Version number of WCSDVARR extension'), 1019 ('D2IM1.NAXES', 2.0, 'Number of independent variables in D2IM function'), 1020 ('D2IM1.AXIS.1', 1.0, 'Axis number of the 1st variable in a D2IM function'), 1021 ('D2IM1.AXIS.2', 2.0, 'Axis number of the 2nd variable in a D2IM function'), 1022 ('D2IMDIS2', 'LOOKUP', 'Detector to image correction type'), 1023 ('D2IM2.EXTVER', 2.0, 'Version number of WCSDVARR extension'), 1024 ('D2IM2.NAXES', 2.0, 'Number of independent variables in D2IM function'), 1025 ('D2IM2.AXIS.1', 1.0, 'Axis number of the 1st variable in a D2IM function'), 1026 ('D2IM2.AXIS.2', 2.0, 'Axis number of the 2nd variable in a D2IM function'), 1027 # ('D2IMERR1', 0.049, 'Maximum error of D2IM correction for axis 1'), 1028 # ('D2IMERR2', 0.035, 'Maximum error of D2IM correction for axis 2'), 1029 # ('D2IMEXT', 'iref$y7b1516hi_d2i.fits', ''), 1030 ] 1031 1032 assert len(wcsdict) == len(refcards) 1033 1034 for k, v, c in refcards: 1035 assert wcsdict[k] == (v, c) 1036 1037 1038def test_sip_broken(): 1039 # This header caused wcslib to segfault because it has a SIP 1040 # specification in a non-default keyword 1041 hdr = get_pkg_data_contents("data/sip-broken.hdr") 1042 1043 wcs.WCS(hdr) 1044 1045 1046def test_no_truncate_crval(): 1047 """ 1048 Regression test for https://github.com/astropy/astropy/issues/4612 1049 """ 1050 w = wcs.WCS(naxis=3) 1051 w.wcs.crval = [50, 50, 2.12345678e11] 1052 w.wcs.cdelt = [1e-3, 1e-3, 1e8] 1053 w.wcs.ctype = ['RA---TAN', 'DEC--TAN', 'FREQ'] 1054 w.wcs.set() 1055 1056 header = w.to_header() 1057 for ii in range(3): 1058 assert header[f'CRVAL{ii + 1}'] == w.wcs.crval[ii] 1059 assert header[f'CDELT{ii + 1}'] == w.wcs.cdelt[ii] 1060 1061 1062def test_no_truncate_crval_try2(): 1063 """ 1064 Regression test for https://github.com/astropy/astropy/issues/4612 1065 """ 1066 w = wcs.WCS(naxis=3) 1067 w.wcs.crval = [50, 50, 2.12345678e11] 1068 w.wcs.cdelt = [1e-5, 1e-5, 1e5] 1069 w.wcs.ctype = ['RA---SIN', 'DEC--SIN', 'FREQ'] 1070 w.wcs.cunit = ['deg', 'deg', 'Hz'] 1071 w.wcs.crpix = [1, 1, 1] 1072 w.wcs.restfrq = 2.34e11 1073 w.wcs.set() 1074 1075 header = w.to_header() 1076 for ii in range(3): 1077 assert header[f'CRVAL{ii + 1}'] == w.wcs.crval[ii] 1078 assert header[f'CDELT{ii + 1}'] == w.wcs.cdelt[ii] 1079 1080 1081def test_no_truncate_crval_p17(): 1082 """ 1083 Regression test for https://github.com/astropy/astropy/issues/5162 1084 """ 1085 w = wcs.WCS(naxis=2) 1086 w.wcs.crval = [50.1234567890123456, 50.1234567890123456] 1087 w.wcs.cdelt = [1e-3, 1e-3] 1088 w.wcs.ctype = ['RA---TAN', 'DEC--TAN'] 1089 w.wcs.set() 1090 1091 header = w.to_header() 1092 assert header['CRVAL1'] != w.wcs.crval[0] 1093 assert header['CRVAL2'] != w.wcs.crval[1] 1094 header = w.to_header(relax=wcs.WCSHDO_P17) 1095 assert header['CRVAL1'] == w.wcs.crval[0] 1096 assert header['CRVAL2'] == w.wcs.crval[1] 1097 1098 1099def test_no_truncate_using_compare(): 1100 """ 1101 Regression test for https://github.com/astropy/astropy/issues/4612 1102 1103 This one uses WCS.wcs.compare and some slightly different values 1104 """ 1105 w = wcs.WCS(naxis=3) 1106 w.wcs.crval = [2.409303333333E+02, 50, 2.12345678e11] 1107 w.wcs.cdelt = [1e-3, 1e-3, 1e8] 1108 w.wcs.ctype = ['RA---TAN', 'DEC--TAN', 'FREQ'] 1109 w.wcs.set() 1110 w2 = wcs.WCS(w.to_header()) 1111 w.wcs.compare(w2.wcs) 1112 1113 1114def test_passing_ImageHDU(): 1115 """ 1116 Passing ImageHDU or PrimaryHDU and comparing it with 1117 wcs initialized from header. For #4493. 1118 """ 1119 path = get_pkg_data_filename('data/validate.fits') 1120 with fits.open(path) as hdulist: 1121 with pytest.warns(wcs.FITSFixedWarning): 1122 wcs_hdu = wcs.WCS(hdulist[0]) 1123 wcs_header = wcs.WCS(hdulist[0].header) 1124 assert wcs_hdu.wcs.compare(wcs_header.wcs) 1125 wcs_hdu = wcs.WCS(hdulist[1]) 1126 wcs_header = wcs.WCS(hdulist[1].header) 1127 assert wcs_hdu.wcs.compare(wcs_header.wcs) 1128 1129 1130def test_inconsistent_sip(): 1131 """ 1132 Test for #4814 1133 """ 1134 hdr = get_pkg_data_contents("data/sip-broken.hdr") 1135 with pytest.warns(None) as wrng: 1136 w = wcs.WCS(hdr) 1137 _check_v71_dateref_warnings(wrng) 1138 with pytest.warns(AstropyWarning): 1139 newhdr = w.to_header(relax=None) 1140 # CTYPE should not include "-SIP" if relax is None 1141 with pytest.warns(None) as wrng: 1142 wnew = wcs.WCS(newhdr) 1143 _check_v71_dateref_warnings(wrng) 1144 assert all(not ctyp.endswith('-SIP') for ctyp in wnew.wcs.ctype) 1145 newhdr = w.to_header(relax=False) 1146 assert('A_0_2' not in newhdr) 1147 # CTYPE should not include "-SIP" if relax is False 1148 with pytest.warns(None) as wrng: 1149 wnew = wcs.WCS(newhdr) 1150 _check_v71_dateref_warnings(wrng) 1151 assert all(not ctyp.endswith('-SIP') for ctyp in wnew.wcs.ctype) 1152 with pytest.warns(AstropyWarning): 1153 newhdr = w.to_header(key="C") 1154 assert('A_0_2' not in newhdr) 1155 # Test writing header with a different key 1156 with pytest.warns(None) as wrng: 1157 wnew = wcs.WCS(newhdr, key='C') 1158 _check_v71_dateref_warnings(wrng) 1159 assert all(not ctyp.endswith('-SIP') for ctyp in wnew.wcs.ctype) 1160 with pytest.warns(AstropyWarning): 1161 newhdr = w.to_header(key=" ") 1162 # Test writing a primary WCS to header 1163 with pytest.warns(None) as wrng: 1164 wnew = wcs.WCS(newhdr) 1165 _check_v71_dateref_warnings(wrng) 1166 assert all(not ctyp.endswith('-SIP') for ctyp in wnew.wcs.ctype) 1167 # Test that "-SIP" is kept into CTYPE if relax=True and 1168 # "-SIP" was in the original header 1169 newhdr = w.to_header(relax=True) 1170 with pytest.warns(None) as wrng: 1171 wnew = wcs.WCS(newhdr) 1172 _check_v71_dateref_warnings(wrng) 1173 assert all(ctyp.endswith('-SIP') for ctyp in wnew.wcs.ctype) 1174 assert('A_0_2' in newhdr) 1175 # Test that SIP coefficients are also written out. 1176 assert wnew.sip is not None 1177 # ######### broken header ########### 1178 # Test that "-SIP" is added to CTYPE if relax=True and 1179 # "-SIP" was not in the original header but SIP coefficients 1180 # are present. 1181 with pytest.warns(None) as wrng: 1182 w = wcs.WCS(hdr) 1183 _check_v71_dateref_warnings(wrng) 1184 w.wcs.ctype = ['RA---TAN', 'DEC--TAN'] 1185 newhdr = w.to_header(relax=True) 1186 with pytest.warns(None) as wrng: 1187 wnew = wcs.WCS(newhdr) 1188 _check_v71_dateref_warnings(wrng) 1189 assert all(ctyp.endswith('-SIP') for ctyp in wnew.wcs.ctype) 1190 1191 1192def test_bounds_check(): 1193 """Test for #4957""" 1194 w = wcs.WCS(naxis=2) 1195 w.wcs.ctype = ["RA---CAR", "DEC--CAR"] 1196 w.wcs.cdelt = [10, 10] 1197 w.wcs.crval = [-90, 90] 1198 w.wcs.crpix = [1, 1] 1199 w.wcs.bounds_check(False, False) 1200 ra, dec = w.wcs_pix2world(300, 0, 0) 1201 assert_allclose(ra, -180) 1202 assert_allclose(dec, -30) 1203 1204 1205def test_naxis(): 1206 w = wcs.WCS(naxis=2) 1207 w.wcs.crval = [1, 1] 1208 w.wcs.cdelt = [0.1, 0.1] 1209 w.wcs.crpix = [1, 1] 1210 w._naxis = [1000, 500] 1211 assert w.pixel_shape == (1000, 500) 1212 assert w.array_shape == (500, 1000) 1213 1214 w.pixel_shape = (99, 59) 1215 assert w._naxis == [99, 59] 1216 1217 w.array_shape = (45, 23) 1218 assert w._naxis == [23, 45] 1219 assert w.pixel_shape == (23, 45) 1220 1221 w.pixel_shape = None 1222 assert w.pixel_bounds is None 1223 1224 1225def test_sip_with_altkey(): 1226 """ 1227 Test that when creating a WCS object using a key, CTYPE with 1228 that key is looked at and not the primary CTYPE. 1229 fix for #5443. 1230 """ 1231 with fits.open(get_pkg_data_filename('data/sip.fits')) as f: 1232 with pytest.warns(wcs.FITSFixedWarning): 1233 w = wcs.WCS(f[0].header) 1234 # create a header with two WCSs. 1235 h1 = w.to_header(relax=True, key='A') 1236 h2 = w.to_header(relax=False) 1237 h1['CTYPE1A'] = "RA---SIN-SIP" 1238 h1['CTYPE2A'] = "DEC--SIN-SIP" 1239 h1.update(h2) 1240 with pytest.warns(None) as wrng: 1241 w = wcs.WCS(h1, key='A') 1242 _check_v71_dateref_warnings(wrng) 1243 assert (w.wcs.ctype == np.array(['RA---SIN-SIP', 'DEC--SIN-SIP'])).all() 1244 1245 1246def test_to_fits_1(): 1247 """ 1248 Test to_fits() with LookupTable distortion. 1249 """ 1250 fits_name = get_pkg_data_filename('data/dist.fits') 1251 with pytest.warns(AstropyDeprecationWarning): 1252 w = wcs.WCS(fits_name) 1253 wfits = w.to_fits() 1254 assert isinstance(wfits, fits.HDUList) 1255 assert isinstance(wfits[0], fits.PrimaryHDU) 1256 assert isinstance(wfits[1], fits.ImageHDU) 1257 1258 1259def test_keyedsip(): 1260 """ 1261 Test sip reading with extra key. 1262 """ 1263 hdr_name = get_pkg_data_filename('data/sip-broken.hdr') 1264 header = fits.Header.fromfile(hdr_name) 1265 del header["CRPIX1"] 1266 del header["CRPIX2"] 1267 1268 w = wcs.WCS(header=header, key="A") 1269 assert isinstance(w.sip, wcs.Sip) 1270 assert w.sip.crpix[0] == 2048 1271 assert w.sip.crpix[1] == 1026 1272 1273 1274def test_zero_size_input(): 1275 with fits.open(get_pkg_data_filename('data/sip.fits')) as f: 1276 with pytest.warns(wcs.FITSFixedWarning): 1277 w = wcs.WCS(f[0].header) 1278 1279 inp = np.zeros((0, 2)) 1280 assert_array_equal(inp, w.all_pix2world(inp, 0)) 1281 assert_array_equal(inp, w.all_world2pix(inp, 0)) 1282 1283 inp = [], [1] 1284 result = w.all_pix2world([], [1], 0) 1285 assert_array_equal(inp[0], result[0]) 1286 assert_array_equal(inp[1], result[1]) 1287 1288 result = w.all_world2pix([], [1], 0) 1289 assert_array_equal(inp[0], result[0]) 1290 assert_array_equal(inp[1], result[1]) 1291 1292 1293def test_scalar_inputs(): 1294 """ 1295 Issue #7845 1296 """ 1297 wcsobj = wcs.WCS(naxis=1) 1298 result = wcsobj.all_pix2world(2, 1) 1299 assert_array_equal(result, [np.array(2.)]) 1300 assert result[0].shape == () 1301 1302 result = wcsobj.all_pix2world([2], 1) 1303 assert_array_equal(result, [np.array([2.])]) 1304 assert result[0].shape == (1,) 1305 1306 1307# Ignore RuntimeWarning raised on s390. 1308@pytest.mark.filterwarnings('ignore:.*invalid value encountered in.*') 1309def test_footprint_contains(): 1310 """ 1311 Test WCS.footprint_contains(skycoord) 1312 """ 1313 1314 header = """ 1315WCSAXES = 2 / Number of coordinate axes 1316CRPIX1 = 1045.0 / Pixel coordinate of reference point 1317CRPIX2 = 1001.0 / Pixel coordinate of reference point 1318PC1_1 = -0.00556448550786 / Coordinate transformation matrix element 1319PC1_2 = -0.001042120133257 / Coordinate transformation matrix element 1320PC2_1 = 0.001181477028705 / Coordinate transformation matrix element 1321PC2_2 = -0.005590809742987 / Coordinate transformation matrix element 1322CDELT1 = 1.0 / [deg] Coordinate increment at reference point 1323CDELT2 = 1.0 / [deg] Coordinate increment at reference point 1324CUNIT1 = 'deg' / Units of coordinate increment and value 1325CUNIT2 = 'deg' / Units of coordinate increment and value 1326CTYPE1 = 'RA---TAN' / TAN (gnomonic) projection + SIP distortions 1327CTYPE2 = 'DEC--TAN' / TAN (gnomonic) projection + SIP distortions 1328CRVAL1 = 250.34971683647 / [deg] Coordinate value at reference point 1329CRVAL2 = 2.2808772582495 / [deg] Coordinate value at reference point 1330LONPOLE = 180.0 / [deg] Native longitude of celestial pole 1331LATPOLE = 2.2808772582495 / [deg] Native latitude of celestial pole 1332RADESYS = 'ICRS' / Equatorial coordinate system 1333MJD-OBS = 58612.339199259 / [d] MJD of observation matching DATE-OBS 1334DATE-OBS= '2019-05-09T08:08:26.816Z' / ISO-8601 observation date matching MJD-OB 1335NAXIS = 2 / NAXIS 1336NAXIS1 = 2136 / length of first array dimension 1337NAXIS2 = 2078 / length of second array dimension 1338 """ # noqa 1339 1340 header = fits.Header.fromstring(header.strip(), '\n') 1341 test_wcs = wcs.WCS(header) 1342 1343 hasCoord = test_wcs.footprint_contains(SkyCoord(254, 2, unit='deg')) 1344 assert hasCoord 1345 1346 hasCoord = test_wcs.footprint_contains(SkyCoord(240, 2, unit='deg')) 1347 assert not hasCoord 1348 1349 hasCoord = test_wcs.footprint_contains(SkyCoord(24, 2, unit='deg')) 1350 assert not hasCoord 1351 1352 1353def test_cunit(): 1354 # Initializing WCS 1355 w1 = wcs.WCS(naxis=2) 1356 w2 = wcs.WCS(naxis=2) 1357 w3 = wcs.WCS(naxis=2) 1358 w4 = wcs.WCS(naxis=2) 1359 # Initializing the values of cunit 1360 w1.wcs.cunit = ['deg', 'm/s'] 1361 w2.wcs.cunit = ['km/h', 'km/h'] 1362 w3.wcs.cunit = ['deg', 'm/s'] 1363 w4.wcs.cunit = ['deg', 'deg'] 1364 1365 # Equality checking a cunit with itself 1366 assert w1.wcs.cunit == w1.wcs.cunit 1367 assert not w1.wcs.cunit != w1.wcs.cunit 1368 # Equality checking of two different cunit object having same values 1369 assert w1.wcs.cunit == w3.wcs.cunit 1370 assert not w1.wcs.cunit != w3.wcs.cunit 1371 # Equality checking of two different cunit object having the same first unit 1372 # but different second unit (see #9154) 1373 assert not w1.wcs.cunit == w4.wcs.cunit 1374 assert w1.wcs.cunit != w4.wcs.cunit 1375 # Inequality checking of two different cunit object having different values 1376 assert not w1.wcs.cunit == w2.wcs.cunit 1377 assert w1.wcs.cunit != w2.wcs.cunit 1378 # Inequality checking of cunit with a list of literals 1379 assert not w1.wcs.cunit == [1, 2, 3] 1380 assert w1.wcs.cunit != [1, 2, 3] 1381 # Inequality checking with some characters 1382 assert not w1.wcs.cunit == ['a', 'b', 'c'] 1383 assert w1.wcs.cunit != ['a', 'b', 'c'] 1384 # Comparison is not implemented TypeError will raise 1385 with pytest.raises(TypeError): 1386 w1.wcs.cunit < w2.wcs.cunit 1387 1388 1389class TestWcsWithTime: 1390 def setup(self): 1391 if _WCSLIB_VER >= Version('7.1'): 1392 fname = get_pkg_data_filename('data/header_with_time_wcslib71.fits') 1393 else: 1394 fname = get_pkg_data_filename('data/header_with_time.fits') 1395 self.header = fits.Header.fromfile(fname) 1396 with pytest.warns(wcs.FITSFixedWarning): 1397 self.w = wcs.WCS(self.header, key='A') 1398 1399 def test_keywods2wcsprm(self): 1400 """ Make sure Wcsprm is populated correctly from the header.""" 1401 1402 ctype = [self.header[val] for val in self.header["CTYPE*"]] 1403 crval = [self.header[val] for val in self.header["CRVAL*"]] 1404 crpix = [self.header[val] for val in self.header["CRPIX*"]] 1405 cdelt = [self.header[val] for val in self.header["CDELT*"]] 1406 cunit = [self.header[val] for val in self.header["CUNIT*"]] 1407 assert list(self.w.wcs.ctype) == ctype 1408 assert list(self.w.wcs.axis_types) == [2200, 2201, 3300, 0] 1409 assert_allclose(self.w.wcs.crval, crval) 1410 assert_allclose(self.w.wcs.crpix, crpix) 1411 assert_allclose(self.w.wcs.cdelt, cdelt) 1412 assert list(self.w.wcs.cunit) == cunit 1413 1414 naxis = self.w.naxis 1415 assert naxis == 4 1416 pc = np.zeros((naxis, naxis), dtype=np.float64) 1417 for i in range(1, 5): 1418 for j in range(1, 5): 1419 if i == j: 1420 pc[i-1, j-1] = self.header.get(f'PC{i}_{j}A', 1) 1421 else: 1422 pc[i-1, j-1] = self.header.get(f'PC{i}_{j}A', 0) 1423 assert_allclose(self.w.wcs.pc, pc) 1424 1425 char_keys = ['timesys', 'trefpos', 'trefdir', 'plephem', 'timeunit', 1426 'dateref', 'dateobs', 'datebeg', 'dateavg', 'dateend'] 1427 for key in char_keys: 1428 assert getattr(self.w.wcs, key) == self.header.get(key, "") 1429 1430 num_keys = ['mjdref', 'mjdobs', 'mjdbeg', 'mjdend', 1431 'jepoch', 'bepoch', 'tstart', 'tstop', 'xposure', 1432 'timsyer', 'timrder', 'timedel', 'timepixr', 1433 'timeoffs', 'telapse', 'czphs', 'cperi'] 1434 1435 for key in num_keys: 1436 if key.upper() == 'MJDREF': 1437 hdrv = [self.header.get('MJDREFIA', np.nan), 1438 self.header.get('MJDREFFA', np.nan)] 1439 else: 1440 hdrv = self.header.get(key, np.nan) 1441 assert_allclose(getattr(self.w.wcs, key), hdrv) 1442 1443 def test_transforms(self): 1444 assert_allclose(self.w.all_pix2world(*self.w.wcs.crpix, 1), 1445 self.w.wcs.crval) 1446 1447 1448def test_invalid_coordinate_masking(): 1449 1450 # Regression test for an issue which caused all coordinates to be set to NaN 1451 # after a transformation rather than just the invalid ones as reported by 1452 # WCSLIB. A specific example of this is that when considering an all-sky 1453 # spectral cube with a spectral axis that is not correlated with the sky 1454 # axes, if transforming pixel coordinates that did not fall 'in' the sky, 1455 # the spectral world value was also masked even though that coordinate 1456 # was valid. 1457 1458 w = wcs.WCS(naxis=3) 1459 w.wcs.ctype = 'VELO_LSR', 'GLON-CAR', 'GLAT-CAR' 1460 w.wcs.crval = -20, 0, 0 1461 w.wcs.crpix = 1, 1441, 241 1462 w.wcs.cdelt = 1.3, -0.125, 0.125 1463 1464 px = [-10, -10, 20] 1465 py = [-10, 10, 20] 1466 pz = [-10, 10, 20] 1467 1468 wx, wy, wz = w.wcs_pix2world(px, py, pz, 0) 1469 1470 # Before fixing this, wx used to return np.nan for the first element 1471 1472 assert_allclose(wx, [-33, -33, 6]) 1473 assert_allclose(wy, [np.nan, 178.75, 177.5]) 1474 assert_allclose(wz, [np.nan, -28.75, -27.5]) 1475 1476 1477def test_no_pixel_area(): 1478 w = wcs.WCS(naxis=3) 1479 1480 # Pixel area cannot be computed 1481 with pytest.raises(ValueError, match='Pixel area is defined only for 2D pixels'): 1482 w.proj_plane_pixel_area() 1483 1484 # Pixel scales still possible 1485 assert_quantity_allclose(w.proj_plane_pixel_scales(), 1) 1486 1487 1488def test_distortion_header(tmpdir): 1489 """ 1490 Test that plate distortion model is correctly described by `wcs.to_header()` 1491 and preserved when creating a Cutout2D from the image, writing it to FITS, 1492 and reading it back from the file. 1493 """ 1494 path = get_pkg_data_filename("data/dss.14.29.56-62.41.05.fits.gz") 1495 cen = np.array((50, 50)) 1496 siz = np.array((20, 20)) 1497 1498 with fits.open(path) as hdulist: 1499 with pytest.warns(wcs.FITSFixedWarning): 1500 w = wcs.WCS(hdulist[0].header) 1501 cut = Cutout2D(hdulist[0].data, position=cen, size=siz, wcs=w) 1502 1503 # This converts the DSS plate solution model with AMD[XY]n coefficients into a 1504 # Template Polynomial Distortion model (TPD.FWD.n coefficients); 1505 # not testing explicitly for the header keywords here. 1506 1507 if _WCSLIB_VER < Version("7.4"): 1508 with pytest.warns(AstropyWarning, match="WCS contains a TPD distortion model in CQDIS"): 1509 w0 = wcs.WCS(w.to_header_string()) 1510 with pytest.warns(AstropyWarning, match="WCS contains a TPD distortion model in CQDIS"): 1511 w1 = wcs.WCS(cut.wcs.to_header_string()) 1512 if _WCSLIB_VER >= Version("7.1"): 1513 pytest.xfail("TPD coefficients incomplete with WCSLIB >= 7.1 < 7.4") 1514 else: 1515 w0 = wcs.WCS(w.to_header_string()) 1516 w1 = wcs.WCS(cut.wcs.to_header_string()) 1517 1518 assert w.pixel_to_world(0, 0).separation(w0.pixel_to_world(0, 0)) < 1.e-3 * u.mas 1519 assert w.pixel_to_world(*cen).separation(w0.pixel_to_world(*cen)) < 1.e-3 * u.mas 1520 1521 assert w.pixel_to_world(*cen).separation(w1.pixel_to_world(*(siz / 2))) < 1.e-3 * u.mas 1522 1523 cutfile = str(tmpdir.join('cutout.fits')) 1524 fits.writeto(cutfile, cut.data, cut.wcs.to_header()) 1525 1526 with fits.open(cutfile) as hdulist: 1527 w2 = wcs.WCS(hdulist[0].header) 1528 1529 assert w.pixel_to_world(*cen).separation(w2.pixel_to_world(*(siz / 2))) < 1.e-3 * u.mas 1530 1531 1532def test_pixlist_wcs_colsel(): 1533 """ 1534 Test selection of a specific pixel list WCS using ``colsel``. See #11412. 1535 """ 1536 hdr_file = get_pkg_data_filename('data/chandra-pixlist-wcs.hdr') 1537 hdr = fits.Header.fromtextfile(hdr_file) 1538 with pytest.warns(wcs.FITSFixedWarning): 1539 w = wcs.WCS(hdr, keysel=['image', 'pixel'], colsel=[11, 12]) 1540 assert w.naxis == 2 1541 assert list(w.wcs.ctype) == ['RA---TAN', 'DEC--TAN'] 1542 assert np.allclose(w.wcs.crval, [229.38051931869, -58.81108068885]) 1543 assert np.allclose(w.wcs.pc, [[1, 0], [0, 1]]) 1544 assert np.allclose(w.wcs.cdelt, [-0.00013666666666666, 0.00013666666666666]) 1545 assert np.allclose(w.wcs.lonpole, 180.) 1546