1# Licensed under a 3-clause BSD style license - see LICENSE.rst 2 3import textwrap 4 5import numpy as np 6import pytest 7 8from astropy.io import fits 9from astropy.nddata.nduncertainty import ( 10 StdDevUncertainty, MissingDataAssociationException, VarianceUncertainty, 11 InverseVariance) 12from astropy import units as u 13from astropy import log 14from astropy.wcs import WCS, FITSFixedWarning 15from astropy.utils import NumpyRNGContext 16from astropy.utils.data import (get_pkg_data_filename, get_pkg_data_filenames, 17 get_pkg_data_contents) 18from astropy.utils.exceptions import AstropyWarning 19 20from astropy.nddata.ccddata import CCDData 21from astropy.nddata import _testing as nd_testing 22from astropy.table import Table 23 24DEFAULT_DATA_SIZE = 100 25 26with NumpyRNGContext(123): 27 _random_array = np.random.normal(size=[DEFAULT_DATA_SIZE, DEFAULT_DATA_SIZE]) 28 29 30def create_ccd_data(): 31 """ 32 Return a CCDData object of size DEFAULT_DATA_SIZE x DEFAULT_DATA_SIZE 33 with units of ADU. 34 """ 35 data = _random_array.copy() 36 fake_meta = {'my_key': 42, 'your_key': 'not 42'} 37 ccd = CCDData(data, unit=u.adu) 38 ccd.header = fake_meta 39 return ccd 40 41 42def test_ccddata_empty(): 43 with pytest.raises(TypeError): 44 CCDData() # empty initializer should fail 45 46 47def test_ccddata_must_have_unit(): 48 with pytest.raises(ValueError): 49 CCDData(np.zeros([2, 2])) 50 51 52def test_ccddata_unit_cannot_be_set_to_none(): 53 ccd_data = create_ccd_data() 54 with pytest.raises(TypeError): 55 ccd_data.unit = None 56 57 58def test_ccddata_meta_header_conflict(): 59 with pytest.raises(ValueError) as exc: 60 CCDData([1, 2, 3], unit='', meta={1: 1}, header={2: 2}) 61 assert "can't have both header and meta." in str(exc.value) 62 63 64def test_ccddata_simple(): 65 ccd_data = create_ccd_data() 66 assert ccd_data.shape == (DEFAULT_DATA_SIZE, DEFAULT_DATA_SIZE) 67 assert ccd_data.size == DEFAULT_DATA_SIZE * DEFAULT_DATA_SIZE 68 assert ccd_data.dtype == np.dtype(float) 69 70 71def test_ccddata_init_with_string_electron_unit(): 72 ccd = CCDData(np.zeros([2, 2]), unit="electron") 73 assert ccd.unit is u.electron 74 75 76def test_initialize_from_FITS(tmpdir): 77 ccd_data = create_ccd_data() 78 hdu = fits.PrimaryHDU(ccd_data) 79 hdulist = fits.HDUList([hdu]) 80 filename = tmpdir.join('afile.fits').strpath 81 hdulist.writeto(filename) 82 cd = CCDData.read(filename, unit=u.electron) 83 assert cd.shape == (DEFAULT_DATA_SIZE, DEFAULT_DATA_SIZE) 84 assert cd.size == DEFAULT_DATA_SIZE * DEFAULT_DATA_SIZE 85 assert np.issubdtype(cd.data.dtype, np.floating) 86 for k, v in hdu.header.items(): 87 assert cd.meta[k] == v 88 89 90def test_initialize_from_fits_with_unit_in_header(tmpdir): 91 fake_img = np.zeros([2, 2]) 92 hdu = fits.PrimaryHDU(fake_img) 93 hdu.header['bunit'] = u.adu.to_string() 94 filename = tmpdir.join('afile.fits').strpath 95 hdu.writeto(filename) 96 ccd = CCDData.read(filename) 97 # ccd should pick up the unit adu from the fits header...did it? 98 assert ccd.unit is u.adu 99 100 # An explicit unit in the read overrides any unit in the FITS file 101 ccd2 = CCDData.read(filename, unit="photon") 102 assert ccd2.unit is u.photon 103 104 105def test_initialize_from_fits_with_ADU_in_header(tmpdir): 106 fake_img = np.zeros([2, 2]) 107 hdu = fits.PrimaryHDU(fake_img) 108 hdu.header['bunit'] = 'ADU' 109 filename = tmpdir.join('afile.fits').strpath 110 hdu.writeto(filename) 111 ccd = CCDData.read(filename) 112 # ccd should pick up the unit adu from the fits header...did it? 113 assert ccd.unit is u.adu 114 115 116def test_initialize_from_fits_with_invalid_unit_in_header(tmpdir): 117 hdu = fits.PrimaryHDU(np.ones((2, 2))) 118 hdu.header['bunit'] = 'definetely-not-a-unit' 119 filename = tmpdir.join('afile.fits').strpath 120 hdu.writeto(filename) 121 with pytest.raises(ValueError): 122 CCDData.read(filename) 123 124 125def test_initialize_from_fits_with_technically_invalid_but_not_really(tmpdir): 126 hdu = fits.PrimaryHDU(np.ones((2, 2))) 127 hdu.header['bunit'] = 'ELECTRONS/S' 128 filename = tmpdir.join('afile.fits').strpath 129 hdu.writeto(filename) 130 ccd = CCDData.read(filename) 131 assert ccd.unit == u.electron/u.s 132 133 134def test_initialize_from_fits_with_data_in_different_extension(tmpdir): 135 fake_img = np.arange(4).reshape(2, 2) 136 hdu1 = fits.PrimaryHDU() 137 hdu2 = fits.ImageHDU(fake_img) 138 hdus = fits.HDUList([hdu1, hdu2]) 139 filename = tmpdir.join('afile.fits').strpath 140 hdus.writeto(filename) 141 ccd = CCDData.read(filename, unit='adu') 142 # ccd should pick up the unit adu from the fits header...did it? 143 np.testing.assert_array_equal(ccd.data, fake_img) 144 # check that the header is the combined header 145 assert hdu2.header + hdu1.header == ccd.header 146 147 148def test_initialize_from_fits_with_extension(tmpdir): 149 fake_img1 = np.zeros([2, 2]) 150 fake_img2 = np.arange(4).reshape(2, 2) 151 hdu0 = fits.PrimaryHDU() 152 hdu1 = fits.ImageHDU(fake_img1) 153 hdu2 = fits.ImageHDU(fake_img2) 154 hdus = fits.HDUList([hdu0, hdu1, hdu2]) 155 filename = tmpdir.join('afile.fits').strpath 156 hdus.writeto(filename) 157 ccd = CCDData.read(filename, hdu=2, unit='adu') 158 # ccd should pick up the unit adu from the fits header...did it? 159 np.testing.assert_array_equal(ccd.data, fake_img2) 160 161 162def test_write_unit_to_hdu(): 163 ccd_data = create_ccd_data() 164 ccd_unit = ccd_data.unit 165 hdulist = ccd_data.to_hdu() 166 assert 'bunit' in hdulist[0].header 167 assert hdulist[0].header['bunit'] == ccd_unit.to_string() 168 169 170def test_initialize_from_FITS_bad_keyword_raises_error(tmpdir): 171 # There are two fits.open keywords that are not permitted in ccdproc: 172 # do_not_scale_image_data and scale_back 173 ccd_data = create_ccd_data() 174 filename = tmpdir.join('test.fits').strpath 175 ccd_data.write(filename) 176 177 with pytest.raises(TypeError): 178 CCDData.read(filename, unit=ccd_data.unit, 179 do_not_scale_image_data=True) 180 with pytest.raises(TypeError): 181 CCDData.read(filename, unit=ccd_data.unit, scale_back=True) 182 183 184def test_ccddata_writer(tmpdir): 185 ccd_data = create_ccd_data() 186 filename = tmpdir.join('test.fits').strpath 187 ccd_data.write(filename) 188 189 ccd_disk = CCDData.read(filename, unit=ccd_data.unit) 190 np.testing.assert_array_equal(ccd_data.data, ccd_disk.data) 191 192 193def test_ccddata_meta_is_case_sensitive(): 194 ccd_data = create_ccd_data() 195 key = 'SoMeKEY' 196 ccd_data.meta[key] = 10 197 assert key.lower() not in ccd_data.meta 198 assert key.upper() not in ccd_data.meta 199 assert key in ccd_data.meta 200 201 202def test_ccddata_meta_is_not_fits_header(): 203 ccd_data = create_ccd_data() 204 ccd_data.meta = {'OBSERVER': 'Edwin Hubble'} 205 assert not isinstance(ccd_data.meta, fits.Header) 206 207 208def test_fromMEF(tmpdir): 209 ccd_data = create_ccd_data() 210 hdu = fits.PrimaryHDU(ccd_data) 211 hdu2 = fits.PrimaryHDU(2 * ccd_data.data) 212 hdulist = fits.HDUList(hdu) 213 hdulist.append(hdu2) 214 filename = tmpdir.join('afile.fits').strpath 215 hdulist.writeto(filename) 216 # by default, we reading from the first extension 217 cd = CCDData.read(filename, unit=u.electron) 218 np.testing.assert_array_equal(cd.data, ccd_data.data) 219 # but reading from the second should work too 220 cd = CCDData.read(filename, hdu=1, unit=u.electron) 221 np.testing.assert_array_equal(cd.data, 2 * ccd_data.data) 222 223 224def test_metafromheader(): 225 hdr = fits.header.Header() 226 hdr['observer'] = 'Edwin Hubble' 227 hdr['exptime'] = '3600' 228 229 d1 = CCDData(np.ones((5, 5)), meta=hdr, unit=u.electron) 230 assert d1.meta['OBSERVER'] == 'Edwin Hubble' 231 assert d1.header['OBSERVER'] == 'Edwin Hubble' 232 233 234def test_metafromdict(): 235 dic = {'OBSERVER': 'Edwin Hubble', 'EXPTIME': 3600} 236 d1 = CCDData(np.ones((5, 5)), meta=dic, unit=u.electron) 237 assert d1.meta['OBSERVER'] == 'Edwin Hubble' 238 239 240def test_header2meta(): 241 hdr = fits.header.Header() 242 hdr['observer'] = 'Edwin Hubble' 243 hdr['exptime'] = '3600' 244 245 d1 = CCDData(np.ones((5, 5)), unit=u.electron) 246 d1.header = hdr 247 assert d1.meta['OBSERVER'] == 'Edwin Hubble' 248 assert d1.header['OBSERVER'] == 'Edwin Hubble' 249 250 251def test_metafromstring_fail(): 252 hdr = 'this is not a valid header' 253 with pytest.raises(TypeError): 254 CCDData(np.ones((5, 5)), meta=hdr, unit=u.adu) 255 256 257def test_setting_bad_uncertainty_raises_error(): 258 ccd_data = create_ccd_data() 259 with pytest.raises(TypeError): 260 # Uncertainty is supposed to be an instance of NDUncertainty 261 ccd_data.uncertainty = 10 262 263 264def test_setting_uncertainty_with_array(): 265 ccd_data = create_ccd_data() 266 ccd_data.uncertainty = None 267 fake_uncertainty = np.sqrt(np.abs(ccd_data.data)) 268 ccd_data.uncertainty = fake_uncertainty.copy() 269 np.testing.assert_array_equal(ccd_data.uncertainty.array, fake_uncertainty) 270 271 272def test_setting_uncertainty_wrong_shape_raises_error(): 273 ccd_data = create_ccd_data() 274 with pytest.raises(ValueError): 275 ccd_data.uncertainty = np.zeros([3, 4]) 276 277 278def test_to_hdu(): 279 ccd_data = create_ccd_data() 280 ccd_data.meta = {'observer': 'Edwin Hubble'} 281 fits_hdulist = ccd_data.to_hdu() 282 assert isinstance(fits_hdulist, fits.HDUList) 283 for k, v in ccd_data.meta.items(): 284 assert fits_hdulist[0].header[k] == v 285 np.testing.assert_array_equal(fits_hdulist[0].data, ccd_data.data) 286 287 288def test_copy(): 289 ccd_data = create_ccd_data() 290 ccd_copy = ccd_data.copy() 291 np.testing.assert_array_equal(ccd_copy.data, ccd_data.data) 292 assert ccd_copy.unit == ccd_data.unit 293 assert ccd_copy.meta == ccd_data.meta 294 295 296@pytest.mark.parametrize('operation,affects_uncertainty', [ 297 ("multiply", True), 298 ("divide", True), 299 ]) 300@pytest.mark.parametrize('operand', [ 301 2.0, 302 2 * u.dimensionless_unscaled, 303 2 * u.photon / u.adu, 304 ]) 305@pytest.mark.parametrize('with_uncertainty', [ 306 True, 307 False]) 308def test_mult_div_overload(operand, with_uncertainty, 309 operation, affects_uncertainty): 310 ccd_data = create_ccd_data() 311 if with_uncertainty: 312 ccd_data.uncertainty = StdDevUncertainty(np.ones_like(ccd_data)) 313 method = getattr(ccd_data, operation) 314 np_method = getattr(np, operation) 315 result = method(operand) 316 assert result is not ccd_data 317 assert isinstance(result, CCDData) 318 assert (result.uncertainty is None or 319 isinstance(result.uncertainty, StdDevUncertainty)) 320 try: 321 op_value = operand.value 322 except AttributeError: 323 op_value = operand 324 325 np.testing.assert_array_equal(result.data, 326 np_method(ccd_data.data, op_value)) 327 if with_uncertainty: 328 if affects_uncertainty: 329 np.testing.assert_array_equal(result.uncertainty.array, 330 np_method(ccd_data.uncertainty.array, 331 op_value)) 332 else: 333 np.testing.assert_array_equal(result.uncertainty.array, 334 ccd_data.uncertainty.array) 335 else: 336 assert result.uncertainty is None 337 338 if isinstance(operand, u.Quantity): 339 # Need the "1 *" below to force arguments to be Quantity to work around 340 # astropy/astropy#2377 341 expected_unit = np_method(1 * ccd_data.unit, 1 * operand.unit).unit 342 assert result.unit == expected_unit 343 else: 344 assert result.unit == ccd_data.unit 345 346 347@pytest.mark.parametrize('operation,affects_uncertainty', [ 348 ("add", False), 349 ("subtract", False), 350 ]) 351@pytest.mark.parametrize('operand,expect_failure', [ 352 (2.0, u.UnitsError), # fail--units don't match image 353 (2 * u.dimensionless_unscaled, u.UnitsError), # same 354 (2 * u.adu, False), 355 ]) 356@pytest.mark.parametrize('with_uncertainty', [ 357 True, 358 False]) 359def test_add_sub_overload(operand, expect_failure, with_uncertainty, 360 operation, affects_uncertainty): 361 ccd_data = create_ccd_data() 362 if with_uncertainty: 363 ccd_data.uncertainty = StdDevUncertainty(np.ones_like(ccd_data)) 364 method = getattr(ccd_data, operation) 365 np_method = getattr(np, operation) 366 if expect_failure: 367 with pytest.raises(expect_failure): 368 result = method(operand) 369 return 370 else: 371 result = method(operand) 372 assert result is not ccd_data 373 assert isinstance(result, CCDData) 374 assert (result.uncertainty is None or 375 isinstance(result.uncertainty, StdDevUncertainty)) 376 try: 377 op_value = operand.value 378 except AttributeError: 379 op_value = operand 380 381 np.testing.assert_array_equal(result.data, 382 np_method(ccd_data.data, op_value)) 383 if with_uncertainty: 384 if affects_uncertainty: 385 np.testing.assert_array_equal(result.uncertainty.array, 386 np_method(ccd_data.uncertainty.array, 387 op_value)) 388 else: 389 np.testing.assert_array_equal(result.uncertainty.array, 390 ccd_data.uncertainty.array) 391 else: 392 assert result.uncertainty is None 393 394 if isinstance(operand, u.Quantity): 395 assert (result.unit == ccd_data.unit and result.unit == operand.unit) 396 else: 397 assert result.unit == ccd_data.unit 398 399 400def test_arithmetic_overload_fails(): 401 ccd_data = create_ccd_data() 402 with pytest.raises(TypeError): 403 ccd_data.multiply("five") 404 405 with pytest.raises(TypeError): 406 ccd_data.divide("five") 407 408 with pytest.raises(TypeError): 409 ccd_data.add("five") 410 411 with pytest.raises(TypeError): 412 ccd_data.subtract("five") 413 414 415def test_arithmetic_no_wcs_compare(): 416 ccd = CCDData(np.ones((10, 10)), unit='') 417 assert ccd.add(ccd, compare_wcs=None).wcs is None 418 assert ccd.subtract(ccd, compare_wcs=None).wcs is None 419 assert ccd.multiply(ccd, compare_wcs=None).wcs is None 420 assert ccd.divide(ccd, compare_wcs=None).wcs is None 421 422 423def test_arithmetic_with_wcs_compare(): 424 def return_true(_, __): 425 return True 426 427 wcs1, wcs2 = nd_testing.create_two_equal_wcs(naxis=2) 428 ccd1 = CCDData(np.ones((10, 10)), unit='', wcs=wcs1) 429 ccd2 = CCDData(np.ones((10, 10)), unit='', wcs=wcs2) 430 nd_testing.assert_wcs_seem_equal( 431 ccd1.add(ccd2, compare_wcs=return_true).wcs, 432 wcs1) 433 nd_testing.assert_wcs_seem_equal( 434 ccd1.subtract(ccd2, compare_wcs=return_true).wcs, 435 wcs1) 436 nd_testing.assert_wcs_seem_equal( 437 ccd1.multiply(ccd2, compare_wcs=return_true).wcs, 438 wcs1) 439 nd_testing.assert_wcs_seem_equal( 440 ccd1.divide(ccd2, compare_wcs=return_true).wcs, 441 wcs1) 442 443 444def test_arithmetic_with_wcs_compare_fail(): 445 def return_false(_, __): 446 return False 447 448 ccd1 = CCDData(np.ones((10, 10)), unit='', wcs=WCS()) 449 ccd2 = CCDData(np.ones((10, 10)), unit='', wcs=WCS()) 450 with pytest.raises(ValueError): 451 ccd1.add(ccd2, compare_wcs=return_false) 452 with pytest.raises(ValueError): 453 ccd1.subtract(ccd2, compare_wcs=return_false) 454 with pytest.raises(ValueError): 455 ccd1.multiply(ccd2, compare_wcs=return_false) 456 with pytest.raises(ValueError): 457 ccd1.divide(ccd2, compare_wcs=return_false) 458 459 460def test_arithmetic_overload_ccddata_operand(): 461 ccd_data = create_ccd_data() 462 ccd_data.uncertainty = StdDevUncertainty(np.ones_like(ccd_data)) 463 operand = ccd_data.copy() 464 result = ccd_data.add(operand) 465 assert len(result.meta) == 0 466 np.testing.assert_array_equal(result.data, 467 2 * ccd_data.data) 468 np.testing.assert_array_almost_equal_nulp( 469 result.uncertainty.array, 470 np.sqrt(2) * ccd_data.uncertainty.array 471 ) 472 473 result = ccd_data.subtract(operand) 474 assert len(result.meta) == 0 475 np.testing.assert_array_equal(result.data, 476 0 * ccd_data.data) 477 np.testing.assert_array_almost_equal_nulp( 478 result.uncertainty.array, 479 np.sqrt(2) * ccd_data.uncertainty.array 480 ) 481 482 result = ccd_data.multiply(operand) 483 assert len(result.meta) == 0 484 np.testing.assert_array_equal(result.data, 485 ccd_data.data ** 2) 486 expected_uncertainty = (np.sqrt(2) * np.abs(ccd_data.data) * 487 ccd_data.uncertainty.array) 488 np.testing.assert_allclose(result.uncertainty.array, 489 expected_uncertainty) 490 491 result = ccd_data.divide(operand) 492 assert len(result.meta) == 0 493 np.testing.assert_array_equal(result.data, 494 np.ones_like(ccd_data.data)) 495 expected_uncertainty = (np.sqrt(2) / np.abs(ccd_data.data) * 496 ccd_data.uncertainty.array) 497 np.testing.assert_allclose(result.uncertainty.array, 498 expected_uncertainty) 499 500 501def test_arithmetic_overload_differing_units(): 502 a = np.array([1, 2, 3]) * u.m 503 b = np.array([1, 2, 3]) * u.cm 504 ccddata = CCDData(a) 505 506 # TODO: Could also be parametrized. 507 res = ccddata.add(b) 508 np.testing.assert_array_almost_equal(res.data, np.add(a, b).value) 509 assert res.unit == np.add(a, b).unit 510 511 res = ccddata.subtract(b) 512 np.testing.assert_array_almost_equal(res.data, np.subtract(a, b).value) 513 assert res.unit == np.subtract(a, b).unit 514 515 res = ccddata.multiply(b) 516 np.testing.assert_array_almost_equal(res.data, np.multiply(a, b).value) 517 assert res.unit == np.multiply(a, b).unit 518 519 res = ccddata.divide(b) 520 np.testing.assert_array_almost_equal(res.data, np.divide(a, b).value) 521 assert res.unit == np.divide(a, b).unit 522 523 524def test_arithmetic_add_with_array(): 525 ccd = CCDData(np.ones((3, 3)), unit='') 526 res = ccd.add(np.arange(3)) 527 np.testing.assert_array_equal(res.data, [[1, 2, 3]] * 3) 528 529 ccd = CCDData(np.ones((3, 3)), unit='adu') 530 with pytest.raises(ValueError): 531 ccd.add(np.arange(3)) 532 533 534def test_arithmetic_subtract_with_array(): 535 ccd = CCDData(np.ones((3, 3)), unit='') 536 res = ccd.subtract(np.arange(3)) 537 np.testing.assert_array_equal(res.data, [[1, 0, -1]] * 3) 538 539 ccd = CCDData(np.ones((3, 3)), unit='adu') 540 with pytest.raises(ValueError): 541 ccd.subtract(np.arange(3)) 542 543 544def test_arithmetic_multiply_with_array(): 545 ccd = CCDData(np.ones((3, 3)) * 3, unit=u.m) 546 res = ccd.multiply(np.ones((3, 3)) * 2) 547 np.testing.assert_array_equal(res.data, [[6, 6, 6]] * 3) 548 assert res.unit == ccd.unit 549 550 551def test_arithmetic_divide_with_array(): 552 ccd = CCDData(np.ones((3, 3)), unit=u.m) 553 res = ccd.divide(np.ones((3, 3)) * 2) 554 np.testing.assert_array_equal(res.data, [[0.5, 0.5, 0.5]] * 3) 555 assert res.unit == ccd.unit 556 557 558def test_history_preserved_if_metadata_is_fits_header(tmpdir): 559 fake_img = np.zeros([2, 2]) 560 hdu = fits.PrimaryHDU(fake_img) 561 hdu.header['history'] = 'one' 562 hdu.header['history'] = 'two' 563 hdu.header['history'] = 'three' 564 assert len(hdu.header['history']) == 3 565 tmp_file = tmpdir.join('temp.fits').strpath 566 hdu.writeto(tmp_file) 567 568 ccd_read = CCDData.read(tmp_file, unit="adu") 569 assert ccd_read.header['history'] == hdu.header['history'] 570 571 572def test_infol_logged_if_unit_in_fits_header(tmpdir): 573 ccd_data = create_ccd_data() 574 tmpfile = tmpdir.join('temp.fits') 575 ccd_data.write(tmpfile.strpath) 576 log.setLevel('INFO') 577 explicit_unit_name = "photon" 578 with log.log_to_list() as log_list: 579 _ = CCDData.read(tmpfile.strpath, unit=explicit_unit_name) 580 assert explicit_unit_name in log_list[0].message 581 582 583def test_wcs_attribute(tmpdir): 584 """ 585 Check that WCS attribute gets added to header, and that if a CCDData 586 object is created from a FITS file with a header, and the WCS attribute 587 is modified, then the CCDData object is turned back into an hdu, the 588 WCS object overwrites the old WCS information in the header. 589 """ 590 ccd_data = create_ccd_data() 591 tmpfile = tmpdir.join('temp.fits') 592 # This wcs example is taken from the astropy.wcs docs. 593 wcs = WCS(naxis=2) 594 wcs.wcs.crpix = np.array(ccd_data.shape) / 2 595 wcs.wcs.cdelt = np.array([-0.066667, 0.066667]) 596 wcs.wcs.crval = [0, -90] 597 wcs.wcs.ctype = ["RA---AIR", "DEC--AIR"] 598 wcs.wcs.set_pv([(2, 1, 45.0)]) 599 ccd_data.header = ccd_data.to_hdu()[0].header 600 ccd_data.header.extend(wcs.to_header(), useblanks=False) 601 ccd_data.write(tmpfile.strpath) 602 603 # Get the header length after it has been extended by the WCS keywords 604 original_header_length = len(ccd_data.header) 605 606 ccd_new = CCDData.read(tmpfile.strpath) 607 # WCS attribute should be set for ccd_new 608 assert ccd_new.wcs is not None 609 # WCS attribute should be equal to wcs above. 610 assert ccd_new.wcs.wcs == wcs.wcs 611 612 # Converting CCDData object with wcs to an hdu shouldn't 613 # create duplicate wcs-related entries in the header. 614 ccd_new_hdu = ccd_new.to_hdu()[0] 615 assert len(ccd_new_hdu.header) == original_header_length 616 617 # Making a CCDData with WCS (but not WCS in the header) should lead to 618 # WCS information in the header when it is converted to an HDU. 619 ccd_wcs_not_in_header = CCDData(ccd_data.data, wcs=wcs, unit="adu") 620 hdu = ccd_wcs_not_in_header.to_hdu()[0] 621 wcs_header = wcs.to_header() 622 for k in wcs_header.keys(): 623 # Skip these keywords if they are in the WCS header because they are 624 # not WCS-specific. 625 if k in ['', 'COMMENT', 'HISTORY']: 626 continue 627 # No keyword from the WCS should be in the header. 628 assert k not in ccd_wcs_not_in_header.header 629 # Every keyword in the WCS should be in the header of the HDU 630 assert hdu.header[k] == wcs_header[k] 631 632 # Now check that if WCS of a CCDData is modified, then the CCDData is 633 # converted to an HDU, the WCS keywords in the header are overwritten 634 # with the appropriate keywords from the header. 635 # 636 # ccd_new has a WCS and WCS keywords in the header, so try modifying 637 # the WCS. 638 ccd_new.wcs.wcs.cdelt *= 2 639 ccd_new_hdu_mod_wcs = ccd_new.to_hdu()[0] 640 assert ccd_new_hdu_mod_wcs.header['CDELT1'] == ccd_new.wcs.wcs.cdelt[0] 641 assert ccd_new_hdu_mod_wcs.header['CDELT2'] == ccd_new.wcs.wcs.cdelt[1] 642 643 644def test_wcs_keywords_removed_from_header(): 645 """ 646 Test, for the file included with the nddata tests, that WCS keywords are 647 properly removed from header. 648 """ 649 from astropy.nddata.ccddata import _KEEP_THESE_KEYWORDS_IN_HEADER 650 keepers = set(_KEEP_THESE_KEYWORDS_IN_HEADER) 651 data_file = get_pkg_data_filename('data/sip-wcs.fits') 652 ccd = CCDData.read(data_file) 653 with pytest.warns(AstropyWarning, 654 match=r'Some non-standard WCS keywords were excluded'): 655 wcs_header = ccd.wcs.to_header() 656 assert not (set(wcs_header) & set(ccd.meta) - keepers) 657 658 # Make sure that exceptions are not raised when trying to remove missing 659 # keywords. o4sp040b0_raw.fits of io.fits is missing keyword 'PC1_1'. 660 data_file1 = get_pkg_data_filename('../../io/fits/tests/data/o4sp040b0_raw.fits') 661 with pytest.warns(FITSFixedWarning, match=r"'unitfix' made the change"): 662 ccd = CCDData.read(data_file1, unit='count') 663 664 665def test_wcs_SIP_coefficient_keywords_removed(): 666 # If SIP polynomials are present, check that no more polynomial 667 # coefficients remain in the header. See #8598 668 669 # The SIP paper is ambiguous as to whether keywords like 670 # A_0_0 can appear in the header for a 2nd order or higher 671 # polynomial. The paper clearly says that the corrections 672 # are only for quadratic or higher order, so A_0_0 and the like 673 # should be zero if they are present, but they apparently can be 674 # there (or at least astrometry.net produces them). 675 676 # astropy WCS does not write those coefficients, so they were 677 # not being removed from the header even though they are WCS-related. 678 679 data_file = get_pkg_data_filename('data/sip-wcs.fits') 680 test_keys = ['A_0_0', 'B_0_1'] 681 682 # Make sure the keywords added to this file for testing are there 683 with fits.open(data_file) as hdu: 684 for key in test_keys: 685 assert key in hdu[0].header 686 687 ccd = CCDData.read(data_file) 688 689 # Now the test...the two keywords above should have been removed. 690 for key in test_keys: 691 assert key not in ccd.header 692 693 694@pytest.mark.filterwarnings('ignore') 695def test_wcs_keyword_removal_for_wcs_test_files(): 696 """ 697 Test, for the WCS test files, that keyword removal works as 698 expected. Those cover a much broader range of WCS types than 699 test_wcs_keywords_removed_from_header. 700 701 Includes regression test for #8597 702 """ 703 from astropy.nddata.ccddata import _generate_wcs_and_update_header 704 from astropy.nddata.ccddata import (_KEEP_THESE_KEYWORDS_IN_HEADER, 705 _CDs, _PCs) 706 707 keepers = set(_KEEP_THESE_KEYWORDS_IN_HEADER) 708 wcs_headers = get_pkg_data_filenames('../../wcs/tests/data', 709 pattern='*.hdr') 710 711 for hdr in wcs_headers: 712 # Skip the files that are expected to be bad... 713 if ('invalid' in hdr or 'nonstandard' in hdr or 'segfault' in hdr or 714 'chandra-pixlist-wcs' in hdr): 715 continue 716 header_string = get_pkg_data_contents(hdr) 717 header = fits.Header.fromstring(header_string) 718 719 wcs = WCS(header_string) 720 header_from_wcs = wcs.to_header(relax=True) 721 722 new_header, new_wcs = _generate_wcs_and_update_header(header) 723 new_wcs_header = new_wcs.to_header(relax=True) 724 725 # Make sure all of the WCS-related keywords generated by astropy 726 # have been removed. 727 assert not (set(new_header) & 728 set(new_wcs_header) - 729 keepers) 730 731 # Check that new_header contains no remaining WCS information. 732 # Specifically, check that 733 # 1. The combination of new_header and new_wcs does not contain 734 # both PCi_j and CDi_j keywords. See #8597. 735 736 # Check for 1 737 final_header = new_header + new_wcs_header 738 final_header_set = set(final_header) 739 740 if _PCs & final_header_set: 741 assert not (_CDs & final_header_set) 742 elif _CDs & final_header_set: 743 assert not (_PCs & final_header_set) 744 745 # Check that the new wcs is the same as the old. 746 for k, v in new_wcs_header.items(): 747 if isinstance(v, str): 748 assert header_from_wcs[k] == v 749 else: 750 np.testing.assert_almost_equal(header_from_wcs[k], v) 751 752 753def test_read_wcs_not_creatable(tmpdir): 754 # The following Header can't be converted to a WCS object. See also #6499. 755 hdr_txt_example_WCS = textwrap.dedent(''' 756 SIMPLE = T / Fits standard 757 BITPIX = 16 / Bits per pixel 758 NAXIS = 2 / Number of axes 759 NAXIS1 = 1104 / Axis length 760 NAXIS2 = 4241 / Axis length 761 CRVAL1 = 164.98110962 / Physical value of the reference pixel X 762 CRVAL2 = 44.34089279 / Physical value of the reference pixel Y 763 CRPIX1 = -34.0 / Reference pixel in X (pixel) 764 CRPIX2 = 2041.0 / Reference pixel in Y (pixel) 765 CDELT1 = 0.10380000 / X Scale projected on detector (#/pix) 766 CDELT2 = 0.10380000 / Y Scale projected on detector (#/pix) 767 CTYPE1 = 'RA---TAN' / Pixel coordinate system 768 CTYPE2 = 'WAVELENGTH' / Pixel coordinate system 769 CUNIT1 = 'degree ' / Units used in both CRVAL1 and CDELT1 770 CUNIT2 = 'nm ' / Units used in both CRVAL2 and CDELT2 771 CD1_1 = 0.20760000 / Pixel Coordinate translation matrix 772 CD1_2 = 0.00000000 / Pixel Coordinate translation matrix 773 CD2_1 = 0.00000000 / Pixel Coordinate translation matrix 774 CD2_2 = 0.10380000 / Pixel Coordinate translation matrix 775 C2YPE1 = 'RA---TAN' / Pixel coordinate system 776 C2YPE2 = 'DEC--TAN' / Pixel coordinate system 777 C2NIT1 = 'degree ' / Units used in both C2VAL1 and C2ELT1 778 C2NIT2 = 'degree ' / Units used in both C2VAL2 and C2ELT2 779 RADECSYS= 'FK5 ' / The equatorial coordinate system 780 ''') 781 hdr = fits.Header.fromstring(hdr_txt_example_WCS, sep='\n') 782 hdul = fits.HDUList([fits.PrimaryHDU(np.ones((4241, 1104)), header=hdr)]) 783 filename = tmpdir.join('afile.fits').strpath 784 hdul.writeto(filename) 785 # The hdr cannot be converted to a WCS object because of an 786 # InconsistentAxisTypesError but it should still open the file 787 ccd = CCDData.read(filename, unit='adu') 788 assert ccd.wcs is None 789 790 791def test_header(): 792 ccd_data = create_ccd_data() 793 a = {'Observer': 'Hubble'} 794 ccd = CCDData(ccd_data, header=a) 795 assert ccd.meta == a 796 797 798def test_wcs_arithmetic(): 799 ccd_data = create_ccd_data() 800 wcs = WCS(naxis=2) 801 ccd_data.wcs = wcs 802 result = ccd_data.multiply(1.0) 803 nd_testing.assert_wcs_seem_equal(result.wcs, wcs) 804 805 806@pytest.mark.parametrize('operation', 807 ['multiply', 'divide', 'add', 'subtract']) 808def test_wcs_arithmetic_ccd(operation): 809 ccd_data = create_ccd_data() 810 ccd_data2 = ccd_data.copy() 811 ccd_data.wcs = WCS(naxis=2) 812 method = getattr(ccd_data, operation) 813 result = method(ccd_data2) 814 nd_testing.assert_wcs_seem_equal(result.wcs, ccd_data.wcs) 815 assert ccd_data2.wcs is None 816 817 818def test_wcs_sip_handling(): 819 """ 820 Check whether the ctypes RA---TAN-SIP and DEC--TAN-SIP survive 821 a roundtrip unchanged. 822 """ 823 data_file = get_pkg_data_filename('data/sip-wcs.fits') 824 825 def check_wcs_ctypes(header): 826 expected_wcs_ctypes = { 827 'CTYPE1': 'RA---TAN-SIP', 828 'CTYPE2': 'DEC--TAN-SIP' 829 } 830 831 return [header[k] == v for k, v in expected_wcs_ctypes.items()] 832 833 ccd_original = CCDData.read(data_file) 834 # After initialization the keywords should be in the WCS, not in the 835 # meta. 836 with fits.open(data_file) as raw: 837 good_ctype = check_wcs_ctypes(raw[0].header) 838 assert all(good_ctype) 839 840 ccd_new = ccd_original.to_hdu() 841 good_ctype = check_wcs_ctypes(ccd_new[0].header) 842 assert all(good_ctype) 843 844 # Try converting to header with wcs_relax=False and 845 # the header should contain the CTYPE keywords without 846 # the -SIP 847 848 ccd_no_relax = ccd_original.to_hdu(wcs_relax=False) 849 good_ctype = check_wcs_ctypes(ccd_no_relax[0].header) 850 assert not any(good_ctype) 851 assert ccd_no_relax[0].header['CTYPE1'] == 'RA---TAN' 852 assert ccd_no_relax[0].header['CTYPE2'] == 'DEC--TAN' 853 854 855@pytest.mark.parametrize('operation', 856 ['multiply', 'divide', 'add', 'subtract']) 857def test_mask_arithmetic_ccd(operation): 858 ccd_data = create_ccd_data() 859 ccd_data2 = ccd_data.copy() 860 ccd_data.mask = (ccd_data.data > 0) 861 method = getattr(ccd_data, operation) 862 result = method(ccd_data2) 863 np.testing.assert_equal(result.mask, ccd_data.mask) 864 865 866def test_write_read_multiextensionfits_mask_default(tmpdir): 867 # Test that if a mask is present the mask is saved and loaded by default. 868 ccd_data = create_ccd_data() 869 ccd_data.mask = ccd_data.data > 10 870 filename = tmpdir.join('afile.fits').strpath 871 ccd_data.write(filename) 872 ccd_after = CCDData.read(filename) 873 assert ccd_after.mask is not None 874 np.testing.assert_array_equal(ccd_data.mask, ccd_after.mask) 875 876 877@pytest.mark.parametrize( 878 'uncertainty_type', 879 [StdDevUncertainty, VarianceUncertainty, InverseVariance]) 880def test_write_read_multiextensionfits_uncertainty_default( 881 tmpdir, uncertainty_type): 882 # Test that if a uncertainty is present it is saved and loaded by default. 883 ccd_data = create_ccd_data() 884 ccd_data.uncertainty = uncertainty_type(ccd_data.data * 10) 885 filename = tmpdir.join('afile.fits').strpath 886 ccd_data.write(filename) 887 ccd_after = CCDData.read(filename) 888 assert ccd_after.uncertainty is not None 889 assert type(ccd_after.uncertainty) is uncertainty_type 890 np.testing.assert_array_equal(ccd_data.uncertainty.array, 891 ccd_after.uncertainty.array) 892 893 894@pytest.mark.parametrize( 895 'uncertainty_type', 896 [StdDevUncertainty, VarianceUncertainty, InverseVariance]) 897def test_write_read_multiextensionfits_uncertainty_different_uncertainty_key( 898 tmpdir, uncertainty_type): 899 # Test that if a uncertainty is present it is saved and loaded by default. 900 ccd_data = create_ccd_data() 901 ccd_data.uncertainty = uncertainty_type(ccd_data.data * 10) 902 filename = tmpdir.join('afile.fits').strpath 903 ccd_data.write(filename, key_uncertainty_type='Blah') 904 ccd_after = CCDData.read(filename, key_uncertainty_type='Blah') 905 assert ccd_after.uncertainty is not None 906 assert type(ccd_after.uncertainty) is uncertainty_type 907 np.testing.assert_array_equal(ccd_data.uncertainty.array, 908 ccd_after.uncertainty.array) 909 910 911def test_write_read_multiextensionfits_not(tmpdir): 912 # Test that writing mask and uncertainty can be disabled 913 ccd_data = create_ccd_data() 914 ccd_data.mask = ccd_data.data > 10 915 ccd_data.uncertainty = StdDevUncertainty(ccd_data.data * 10) 916 filename = tmpdir.join('afile.fits').strpath 917 ccd_data.write(filename, hdu_mask=None, hdu_uncertainty=None) 918 ccd_after = CCDData.read(filename) 919 assert ccd_after.uncertainty is None 920 assert ccd_after.mask is None 921 922 923def test_write_read_multiextensionfits_custom_ext_names(tmpdir): 924 # Test writing mask, uncertainty in another extension than default 925 ccd_data = create_ccd_data() 926 ccd_data.mask = ccd_data.data > 10 927 ccd_data.uncertainty = StdDevUncertainty(ccd_data.data * 10) 928 filename = tmpdir.join('afile.fits').strpath 929 ccd_data.write(filename, hdu_mask='Fun', hdu_uncertainty='NoFun') 930 931 # Try reading with defaults extension names 932 ccd_after = CCDData.read(filename) 933 assert ccd_after.uncertainty is None 934 assert ccd_after.mask is None 935 936 # Try reading with custom extension names 937 ccd_after = CCDData.read(filename, hdu_mask='Fun', hdu_uncertainty='NoFun') 938 assert ccd_after.uncertainty is not None 939 assert ccd_after.mask is not None 940 np.testing.assert_array_equal(ccd_data.mask, ccd_after.mask) 941 np.testing.assert_array_equal(ccd_data.uncertainty.array, 942 ccd_after.uncertainty.array) 943 944 945def test_read_old_style_multiextensionfits(tmpdir): 946 # Regression test for https://github.com/astropy/ccdproc/issues/664 947 # 948 # Prior to astropy 3.1 there was no uncertainty type saved 949 # in the multiextension fits files generated by CCDData 950 # because the uncertainty had to be StandardDevUncertainty. 951 # 952 # Current version should be able to read those in. 953 # 954 size = 4 955 # Value of the variables below are not important to the test. 956 data = np.zeros([size, size]) 957 mask = data > 0.9 958 uncert = np.sqrt(data) 959 960 ccd = CCDData(data=data, mask=mask, uncertainty=uncert, unit='adu') 961 # We'll create the file manually to ensure we have the 962 # right extension names and no uncertainty type. 963 hdulist = ccd.to_hdu() 964 del hdulist[2].header['UTYPE'] 965 file_name = tmpdir.join('old_ccddata_mef.fits').strpath 966 hdulist.writeto(file_name) 967 968 ccd = CCDData.read(file_name) 969 970 assert isinstance(ccd.uncertainty, StdDevUncertainty) 971 972 973def test_wcs(): 974 ccd_data = create_ccd_data() 975 wcs = WCS(naxis=2) 976 ccd_data.wcs = wcs 977 assert ccd_data.wcs is wcs 978 979 980def test_recognized_fits_formats_for_read_write(tmpdir): 981 # These are the extensions that are supposed to be supported. 982 ccd_data = create_ccd_data() 983 supported_extensions = ['fit', 'fits', 'fts'] 984 985 for ext in supported_extensions: 986 path = tmpdir.join(f"test.{ext}") 987 ccd_data.write(path.strpath) 988 from_disk = CCDData.read(path.strpath) 989 assert (ccd_data.data == from_disk.data).all() 990 991 992def test_stddevuncertainty_compat_descriptor_no_parent(): 993 with pytest.raises(MissingDataAssociationException): 994 StdDevUncertainty(np.ones((10, 10))).parent_nddata 995 996 997def test_stddevuncertainty_compat_descriptor_no_weakref(): 998 # TODO: Remove this test if astropy 1.0 isn't supported anymore 999 # This test might create a Memoryleak on purpose, so the last lines after 1000 # the assert are IMPORTANT cleanup. 1001 ccd = CCDData(np.ones((10, 10)), unit='') 1002 uncert = StdDevUncertainty(np.ones((10, 10))) 1003 uncert._parent_nddata = ccd 1004 assert uncert.parent_nddata is ccd 1005 uncert._parent_nddata = None 1006 1007 1008# https://github.com/astropy/astropy/issues/7595 1009def test_read_returns_image(tmpdir): 1010 # Test if CCData.read returns a image when reading a fits file containing 1011 # a table and image, in that order. 1012 tbl = Table(np.ones(10).reshape(5, 2)) 1013 img = np.ones((5, 5)) 1014 hdul = fits.HDUList(hdus=[fits.PrimaryHDU(), fits.TableHDU(tbl.as_array()), 1015 fits.ImageHDU(img)]) 1016 filename = tmpdir.join('table_image.fits').strpath 1017 hdul.writeto(filename) 1018 ccd = CCDData.read(filename, unit='adu') 1019 # Expecting to get (5, 5), the size of the image 1020 assert ccd.data.shape == (5, 5) 1021 1022 1023# https://github.com/astropy/astropy/issues/9664 1024def test_sliced_ccdata_to_hdu(): 1025 wcs = WCS(naxis=2) 1026 wcs.wcs.crpix = 10, 10 1027 ccd = CCDData(np.ones((10, 10)), wcs=wcs, unit='pixel') 1028 trimmed = ccd[2:-2, 2:-2] 1029 hdul = trimmed.to_hdu() 1030 assert isinstance(hdul, fits.HDUList) 1031 assert hdul[0].header['CRPIX1'] == 8 1032 assert hdul[0].header['CRPIX2'] == 8 1033