1# Licensed under a 3-clause BSD style license - see LICENSE.rst 2 3import pytest 4 5from io import StringIO 6from itertools import product 7 8from packaging.version import Version 9import numpy as np 10from numpy.testing import assert_almost_equal, assert_equal, assert_allclose 11 12from astropy.utils.data import get_pkg_data_contents, get_pkg_data_filename 13from astropy.utils.exceptions import AstropyUserWarning 14from astropy.time import Time 15from astropy import units as u 16from astropy.utils import unbroadcast 17from astropy.coordinates import SkyCoord, EarthLocation, ITRS 18from astropy.units import Quantity 19from astropy.io import fits 20 21from astropy.wcs import _wcs # noqa 22from astropy.wcs.wcs import (WCS, Sip, WCSSUB_LONGITUDE, WCSSUB_LATITUDE, 23 FITSFixedWarning) 24from astropy.wcs.wcsapi.fitswcs import SlicedFITSWCS 25from astropy.wcs.utils import (proj_plane_pixel_scales, 26 is_proj_plane_distorted, 27 non_celestial_pixel_scales, 28 wcs_to_celestial_frame, 29 celestial_frame_to_wcs, skycoord_to_pixel, 30 pixel_to_skycoord, custom_wcs_to_frame_mappings, 31 custom_frame_to_wcs_mappings, 32 add_stokes_axis_to_wcs, 33 pixel_to_pixel, 34 _split_matrix, 35 _pixel_to_pixel_correlation_matrix, 36 _pixel_to_world_correlation_matrix, 37 local_partial_pixel_derivatives, 38 fit_wcs_from_points, 39 obsgeo_to_frame) 40from astropy.utils.compat.optional_deps import HAS_SCIPY # noqa 41 42 43def test_wcs_dropping(): 44 wcs = WCS(naxis=4) 45 wcs.wcs.pc = np.zeros([4, 4]) 46 np.fill_diagonal(wcs.wcs.pc, np.arange(1, 5)) 47 pc = wcs.wcs.pc # for later use below 48 49 dropped = wcs.dropaxis(0) 50 assert np.all(dropped.wcs.get_pc().diagonal() == np.array([2, 3, 4])) 51 dropped = wcs.dropaxis(1) 52 assert np.all(dropped.wcs.get_pc().diagonal() == np.array([1, 3, 4])) 53 dropped = wcs.dropaxis(2) 54 assert np.all(dropped.wcs.get_pc().diagonal() == np.array([1, 2, 4])) 55 dropped = wcs.dropaxis(3) 56 assert np.all(dropped.wcs.get_pc().diagonal() == np.array([1, 2, 3])) 57 58 wcs = WCS(naxis=4) 59 wcs.wcs.cd = pc 60 61 dropped = wcs.dropaxis(0) 62 assert np.all(dropped.wcs.get_pc().diagonal() == np.array([2, 3, 4])) 63 dropped = wcs.dropaxis(1) 64 assert np.all(dropped.wcs.get_pc().diagonal() == np.array([1, 3, 4])) 65 dropped = wcs.dropaxis(2) 66 assert np.all(dropped.wcs.get_pc().diagonal() == np.array([1, 2, 4])) 67 dropped = wcs.dropaxis(3) 68 assert np.all(dropped.wcs.get_pc().diagonal() == np.array([1, 2, 3])) 69 70 71def test_wcs_swapping(): 72 wcs = WCS(naxis=4) 73 wcs.wcs.pc = np.zeros([4, 4]) 74 np.fill_diagonal(wcs.wcs.pc, np.arange(1, 5)) 75 pc = wcs.wcs.pc # for later use below 76 77 swapped = wcs.swapaxes(0, 1) 78 assert np.all(swapped.wcs.get_pc().diagonal() == np.array([2, 1, 3, 4])) 79 swapped = wcs.swapaxes(0, 3) 80 assert np.all(swapped.wcs.get_pc().diagonal() == np.array([4, 2, 3, 1])) 81 swapped = wcs.swapaxes(2, 3) 82 assert np.all(swapped.wcs.get_pc().diagonal() == np.array([1, 2, 4, 3])) 83 84 wcs = WCS(naxis=4) 85 wcs.wcs.cd = pc 86 87 swapped = wcs.swapaxes(0, 1) 88 assert np.all(swapped.wcs.get_pc().diagonal() == np.array([2, 1, 3, 4])) 89 swapped = wcs.swapaxes(0, 3) 90 assert np.all(swapped.wcs.get_pc().diagonal() == np.array([4, 2, 3, 1])) 91 swapped = wcs.swapaxes(2, 3) 92 assert np.all(swapped.wcs.get_pc().diagonal() == np.array([1, 2, 4, 3])) 93 94 95@pytest.mark.parametrize('ndim', (2, 3)) 96def test_add_stokes(ndim): 97 wcs = WCS(naxis=ndim) 98 99 for ii in range(ndim + 1): 100 outwcs = add_stokes_axis_to_wcs(wcs, ii) 101 assert outwcs.wcs.naxis == ndim + 1 102 assert outwcs.wcs.ctype[ii] == 'STOKES' 103 assert outwcs.wcs.cname[ii] == 'STOKES' 104 105 106def test_slice(): 107 mywcs = WCS(naxis=2) 108 mywcs.wcs.crval = [1, 1] 109 mywcs.wcs.cdelt = [0.1, 0.1] 110 mywcs.wcs.crpix = [1, 1] 111 mywcs._naxis = [1000, 500] 112 pscale = 0.1 # from cdelt 113 114 slice_wcs = mywcs.slice([slice(1, None), slice(0, None)]) 115 assert np.all(slice_wcs.wcs.crpix == np.array([1, 0])) 116 assert slice_wcs._naxis == [1000, 499] 117 118 # test that CRPIX maps to CRVAL: 119 assert_allclose( 120 slice_wcs.wcs_pix2world(*slice_wcs.wcs.crpix, 1), 121 slice_wcs.wcs.crval, rtol=0.0, atol=1e-6 * pscale 122 ) 123 124 slice_wcs = mywcs.slice([slice(1, None, 2), slice(0, None, 4)]) 125 assert np.all(slice_wcs.wcs.crpix == np.array([0.625, 0.25])) 126 assert np.all(slice_wcs.wcs.cdelt == np.array([0.4, 0.2])) 127 assert slice_wcs._naxis == [250, 250] 128 129 slice_wcs = mywcs.slice([slice(None, None, 2), slice(0, None, 2)]) 130 assert np.all(slice_wcs.wcs.cdelt == np.array([0.2, 0.2])) 131 assert slice_wcs._naxis == [500, 250] 132 133 # Non-integral values do not alter the naxis attribute 134 with pytest.warns(AstropyUserWarning): 135 slice_wcs = mywcs.slice([slice(50.), slice(20.)]) 136 assert slice_wcs._naxis == [1000, 500] 137 with pytest.warns(AstropyUserWarning): 138 slice_wcs = mywcs.slice([slice(50.), slice(20)]) 139 assert slice_wcs._naxis == [20, 500] 140 with pytest.warns(AstropyUserWarning): 141 slice_wcs = mywcs.slice([slice(50), slice(20.5)]) 142 assert slice_wcs._naxis == [1000, 50] 143 144 145def test_slice_with_sip(): 146 mywcs = WCS(naxis=2) 147 mywcs.wcs.crval = [1, 1] 148 mywcs.wcs.cdelt = [0.1, 0.1] 149 mywcs.wcs.crpix = [1, 1] 150 mywcs._naxis = [1000, 500] 151 mywcs.wcs.ctype = ['RA---TAN-SIP', 'DEC--TAN-SIP'] 152 a = np.array( 153 [[0, 0, 5.33092692e-08, 3.73753773e-11, -2.02111473e-13], 154 [0, 2.44084308e-05, 2.81394789e-11, 5.17856895e-13, 0.0], 155 [-2.41334657e-07, 1.29289255e-10, 2.35753629e-14, 0.0, 0.0], 156 [-2.37162007e-10, 5.43714947e-13, 0.0, 0.0, 0.0], 157 [-2.81029767e-13, 0.0, 0.0, 0.0, 0.0]] 158 ) 159 b = np.array( 160 [[0, 0, 2.99270374e-05, -2.38136074e-10, 7.23205168e-13], 161 [0, -1.71073858e-07, 6.31243431e-11, -5.16744347e-14, 0.0], 162 [6.95458963e-06, -3.08278961e-10, -1.75800917e-13, 0.0, 0.0], 163 [3.51974159e-11, 5.60993016e-14, 0.0, 0.0, 0.0], 164 [-5.92438525e-13, 0.0, 0.0, 0.0, 0.0]] 165 ) 166 mywcs.sip = Sip(a, b, None, None, mywcs.wcs.crpix) 167 mywcs.wcs.set() 168 pscale = 0.1 # from cdelt 169 170 slice_wcs = mywcs.slice([slice(1, None), slice(0, None)]) 171 # test that CRPIX maps to CRVAL: 172 assert_allclose( 173 slice_wcs.all_pix2world(*slice_wcs.wcs.crpix, 1), 174 slice_wcs.wcs.crval, rtol=0.0, atol=1e-6 * pscale 175 ) 176 177 slice_wcs = mywcs.slice([slice(1, None, 2), slice(0, None, 4)]) 178 # test that CRPIX maps to CRVAL: 179 assert_allclose( 180 slice_wcs.all_pix2world(*slice_wcs.wcs.crpix, 1), 181 slice_wcs.wcs.crval, rtol=0.0, atol=1e-6 * pscale 182 ) 183 184 185def test_slice_getitem(): 186 mywcs = WCS(naxis=2) 187 mywcs.wcs.crval = [1, 1] 188 mywcs.wcs.cdelt = [0.1, 0.1] 189 mywcs.wcs.crpix = [1, 1] 190 191 slice_wcs = mywcs[1::2, 0::4] 192 assert np.all(slice_wcs.wcs.crpix == np.array([0.625, 0.25])) 193 assert np.all(slice_wcs.wcs.cdelt == np.array([0.4, 0.2])) 194 195 mywcs.wcs.crpix = [2, 2] 196 slice_wcs = mywcs[1::2, 0::4] 197 assert np.all(slice_wcs.wcs.crpix == np.array([0.875, 0.75])) 198 assert np.all(slice_wcs.wcs.cdelt == np.array([0.4, 0.2])) 199 200 # Default: numpy order 201 slice_wcs = mywcs[1::2] 202 assert np.all(slice_wcs.wcs.crpix == np.array([2, 0.75])) 203 assert np.all(slice_wcs.wcs.cdelt == np.array([0.1, 0.2])) 204 205 206def test_slice_fitsorder(): 207 mywcs = WCS(naxis=2) 208 mywcs.wcs.crval = [1, 1] 209 mywcs.wcs.cdelt = [0.1, 0.1] 210 mywcs.wcs.crpix = [1, 1] 211 212 slice_wcs = mywcs.slice([slice(1, None), slice(0, None)], numpy_order=False) 213 assert np.all(slice_wcs.wcs.crpix == np.array([0, 1])) 214 215 slice_wcs = mywcs.slice([slice(1, None, 2), slice(0, None, 4)], numpy_order=False) 216 assert np.all(slice_wcs.wcs.crpix == np.array([0.25, 0.625])) 217 assert np.all(slice_wcs.wcs.cdelt == np.array([0.2, 0.4])) 218 219 slice_wcs = mywcs.slice([slice(1, None, 2)], numpy_order=False) 220 assert np.all(slice_wcs.wcs.crpix == np.array([0.25, 1])) 221 assert np.all(slice_wcs.wcs.cdelt == np.array([0.2, 0.1])) 222 223 224def test_slice_wcs(): 225 mywcs = WCS(naxis=2) 226 227 sub = mywcs[0] 228 assert isinstance(sub, SlicedFITSWCS) 229 230 with pytest.raises(IndexError) as exc: 231 mywcs[0, ::2] 232 assert exc.value.args[0] == "Slicing WCS with a step is not supported." 233 234 235def test_axis_names(): 236 mywcs = WCS(naxis=4) 237 mywcs.wcs.ctype = ['RA---TAN', 'DEC--TAN', 'VOPT-LSR', 'STOKES'] 238 239 assert mywcs.axis_type_names == ['RA', 'DEC', 'VOPT', 'STOKES'] 240 241 mywcs.wcs.cname = ['RA', 'DEC', 'VOPT', 'STOKES'] 242 243 assert mywcs.axis_type_names == ['RA', 'DEC', 'VOPT', 'STOKES'] 244 245 246def test_celestial(): 247 mywcs = WCS(naxis=4) 248 mywcs.wcs.ctype = ['RA---TAN', 'DEC--TAN', 'VOPT', 'STOKES'] 249 cel = mywcs.celestial 250 assert tuple(cel.wcs.ctype) == ('RA---TAN', 'DEC--TAN') 251 assert cel.axis_type_names == ['RA', 'DEC'] 252 253 254def test_wcs_to_celestial_frame(): 255 256 # Import astropy.coordinates here to avoid circular imports 257 from astropy.coordinates.builtin_frames import ICRS, ITRS, FK5, FK4, Galactic 258 259 mywcs = WCS(naxis=2) 260 mywcs.wcs.set() 261 with pytest.raises(ValueError, match="Could not determine celestial frame " 262 "corresponding to the specified WCS object"): 263 assert wcs_to_celestial_frame(mywcs) is None 264 265 mywcs = WCS(naxis=2) 266 mywcs.wcs.ctype = ['XOFFSET', 'YOFFSET'] 267 mywcs.wcs.set() 268 with pytest.raises(ValueError): 269 assert wcs_to_celestial_frame(mywcs) is None 270 271 mywcs = WCS(naxis=2) 272 mywcs.wcs.ctype = ['RA---TAN', 'DEC--TAN'] 273 mywcs.wcs.set() 274 frame = wcs_to_celestial_frame(mywcs) 275 assert isinstance(frame, ICRS) 276 277 mywcs = WCS(naxis=2) 278 mywcs.wcs.ctype = ['RA---TAN', 'DEC--TAN'] 279 mywcs.wcs.equinox = 1987. 280 mywcs.wcs.set() 281 print(mywcs.to_header()) 282 frame = wcs_to_celestial_frame(mywcs) 283 assert isinstance(frame, FK5) 284 assert frame.equinox == Time(1987., format='jyear') 285 286 mywcs = WCS(naxis=2) 287 mywcs.wcs.ctype = ['RA---TAN', 'DEC--TAN'] 288 mywcs.wcs.equinox = 1982 289 mywcs.wcs.set() 290 frame = wcs_to_celestial_frame(mywcs) 291 assert isinstance(frame, FK4) 292 assert frame.equinox == Time(1982., format='byear') 293 294 mywcs = WCS(naxis=2) 295 mywcs.wcs.ctype = ['GLON-SIN', 'GLAT-SIN'] 296 mywcs.wcs.set() 297 frame = wcs_to_celestial_frame(mywcs) 298 assert isinstance(frame, Galactic) 299 300 mywcs = WCS(naxis=2) 301 mywcs.wcs.ctype = ['TLON-CAR', 'TLAT-CAR'] 302 mywcs.wcs.dateobs = '2017-08-17T12:41:04.430' 303 mywcs.wcs.set() 304 frame = wcs_to_celestial_frame(mywcs) 305 assert isinstance(frame, ITRS) 306 assert frame.obstime == Time('2017-08-17T12:41:04.430') 307 308 for equinox in [np.nan, 1987, 1982]: 309 mywcs = WCS(naxis=2) 310 mywcs.wcs.ctype = ['RA---TAN', 'DEC--TAN'] 311 mywcs.wcs.radesys = 'ICRS' 312 mywcs.wcs.equinox = equinox 313 mywcs.wcs.set() 314 frame = wcs_to_celestial_frame(mywcs) 315 assert isinstance(frame, ICRS) 316 317 # Flipped order 318 mywcs = WCS(naxis=2) 319 mywcs.wcs.ctype = ['DEC--TAN', 'RA---TAN'] 320 mywcs.wcs.set() 321 frame = wcs_to_celestial_frame(mywcs) 322 assert isinstance(frame, ICRS) 323 324 # More than two dimensions 325 mywcs = WCS(naxis=3) 326 mywcs.wcs.ctype = ['DEC--TAN', 'VELOCITY', 'RA---TAN'] 327 mywcs.wcs.set() 328 frame = wcs_to_celestial_frame(mywcs) 329 assert isinstance(frame, ICRS) 330 331 mywcs = WCS(naxis=3) 332 mywcs.wcs.ctype = ['GLAT-CAR', 'VELOCITY', 'GLON-CAR'] 333 mywcs.wcs.set() 334 frame = wcs_to_celestial_frame(mywcs) 335 assert isinstance(frame, Galactic) 336 337 338def test_wcs_to_celestial_frame_correlated(): 339 340 # Regression test for a bug that caused wcs_to_celestial_frame to fail when 341 # the celestial axes were correlated with other axes. 342 343 # Import astropy.coordinates here to avoid circular imports 344 from astropy.coordinates.builtin_frames import ICRS 345 346 mywcs = WCS(naxis=3) 347 mywcs.wcs.ctype = 'RA---TAN', 'DEC--TAN', 'FREQ' 348 mywcs.wcs.cd = np.ones((3, 3)) 349 mywcs.wcs.set() 350 frame = wcs_to_celestial_frame(mywcs) 351 assert isinstance(frame, ICRS) 352 353 354def test_wcs_to_celestial_frame_extend(): 355 356 mywcs = WCS(naxis=2) 357 mywcs.wcs.ctype = ['XOFFSET', 'YOFFSET'] 358 mywcs.wcs.set() 359 with pytest.raises(ValueError): 360 wcs_to_celestial_frame(mywcs) 361 362 class OffsetFrame: 363 pass 364 365 def identify_offset(wcs): 366 if wcs.wcs.ctype[0].endswith('OFFSET') and wcs.wcs.ctype[1].endswith('OFFSET'): 367 return OffsetFrame() 368 369 with custom_wcs_to_frame_mappings(identify_offset): 370 frame = wcs_to_celestial_frame(mywcs) 371 assert isinstance(frame, OffsetFrame) 372 373 # Check that things are back to normal after the context manager 374 with pytest.raises(ValueError): 375 wcs_to_celestial_frame(mywcs) 376 377 378def test_celestial_frame_to_wcs(): 379 380 # Import astropy.coordinates here to avoid circular imports 381 from astropy.coordinates import ICRS, ITRS, FK5, FK4, FK4NoETerms, Galactic, BaseCoordinateFrame 382 383 class FakeFrame(BaseCoordinateFrame): 384 pass 385 386 frame = FakeFrame() 387 with pytest.raises(ValueError) as exc: 388 celestial_frame_to_wcs(frame) 389 assert exc.value.args[0] == ("Could not determine WCS corresponding to " 390 "the specified coordinate frame.") 391 392 frame = ICRS() 393 mywcs = celestial_frame_to_wcs(frame) 394 mywcs.wcs.set() 395 assert tuple(mywcs.wcs.ctype) == ('RA---TAN', 'DEC--TAN') 396 assert mywcs.wcs.radesys == 'ICRS' 397 assert np.isnan(mywcs.wcs.equinox) 398 assert mywcs.wcs.lonpole == 180 399 assert mywcs.wcs.latpole == 0 400 401 frame = FK5(equinox='J1987') 402 mywcs = celestial_frame_to_wcs(frame) 403 assert tuple(mywcs.wcs.ctype) == ('RA---TAN', 'DEC--TAN') 404 assert mywcs.wcs.radesys == 'FK5' 405 assert mywcs.wcs.equinox == 1987. 406 407 frame = FK4(equinox='B1982') 408 mywcs = celestial_frame_to_wcs(frame) 409 assert tuple(mywcs.wcs.ctype) == ('RA---TAN', 'DEC--TAN') 410 assert mywcs.wcs.radesys == 'FK4' 411 assert mywcs.wcs.equinox == 1982. 412 413 frame = FK4NoETerms(equinox='B1982') 414 mywcs = celestial_frame_to_wcs(frame) 415 assert tuple(mywcs.wcs.ctype) == ('RA---TAN', 'DEC--TAN') 416 assert mywcs.wcs.radesys == 'FK4-NO-E' 417 assert mywcs.wcs.equinox == 1982. 418 419 frame = Galactic() 420 mywcs = celestial_frame_to_wcs(frame) 421 assert tuple(mywcs.wcs.ctype) == ('GLON-TAN', 'GLAT-TAN') 422 assert mywcs.wcs.radesys == '' 423 assert np.isnan(mywcs.wcs.equinox) 424 425 frame = Galactic() 426 mywcs = celestial_frame_to_wcs(frame, projection='CAR') 427 assert tuple(mywcs.wcs.ctype) == ('GLON-CAR', 'GLAT-CAR') 428 assert mywcs.wcs.radesys == '' 429 assert np.isnan(mywcs.wcs.equinox) 430 431 frame = Galactic() 432 mywcs = celestial_frame_to_wcs(frame, projection='CAR') 433 mywcs.wcs.crval = [100, -30] 434 mywcs.wcs.set() 435 assert_allclose((mywcs.wcs.lonpole, mywcs.wcs.latpole), (180, 60)) 436 437 frame = ITRS(obstime=Time('2017-08-17T12:41:04.43')) 438 mywcs = celestial_frame_to_wcs(frame, projection='CAR') 439 assert tuple(mywcs.wcs.ctype) == ('TLON-CAR', 'TLAT-CAR') 440 assert mywcs.wcs.radesys == 'ITRS' 441 assert mywcs.wcs.dateobs == '2017-08-17T12:41:04.430' 442 443 frame = ITRS() 444 mywcs = celestial_frame_to_wcs(frame, projection='CAR') 445 assert tuple(mywcs.wcs.ctype) == ('TLON-CAR', 'TLAT-CAR') 446 assert mywcs.wcs.radesys == 'ITRS' 447 assert mywcs.wcs.dateobs == Time('J2000').utc.fits 448 449 450def test_celestial_frame_to_wcs_extend(): 451 452 class OffsetFrame: 453 pass 454 455 frame = OffsetFrame() 456 457 with pytest.raises(ValueError): 458 celestial_frame_to_wcs(frame) 459 460 def identify_offset(frame, projection=None): 461 if isinstance(frame, OffsetFrame): 462 wcs = WCS(naxis=2) 463 wcs.wcs.ctype = ['XOFFSET', 'YOFFSET'] 464 return wcs 465 466 with custom_frame_to_wcs_mappings(identify_offset): 467 mywcs = celestial_frame_to_wcs(frame) 468 assert tuple(mywcs.wcs.ctype) == ('XOFFSET', 'YOFFSET') 469 470 # Check that things are back to normal after the context manager 471 with pytest.raises(ValueError): 472 celestial_frame_to_wcs(frame) 473 474 475def test_pixscale_nodrop(): 476 mywcs = WCS(naxis=2) 477 mywcs.wcs.cdelt = [0.1, 0.2] 478 mywcs.wcs.ctype = ['RA---TAN', 'DEC--TAN'] 479 assert_almost_equal(proj_plane_pixel_scales(mywcs), (0.1, 0.2)) 480 481 mywcs.wcs.cdelt = [-0.1, 0.2] 482 assert_almost_equal(proj_plane_pixel_scales(mywcs), (0.1, 0.2)) 483 484 485def test_pixscale_withdrop(): 486 mywcs = WCS(naxis=3) 487 mywcs.wcs.cdelt = [0.1, 0.2, 1] 488 mywcs.wcs.ctype = ['RA---TAN', 'DEC--TAN', 'VOPT'] 489 assert_almost_equal(proj_plane_pixel_scales(mywcs.celestial), (0.1, 0.2)) 490 491 mywcs.wcs.cdelt = [-0.1, 0.2, 1] 492 assert_almost_equal(proj_plane_pixel_scales(mywcs.celestial), (0.1, 0.2)) 493 494 495def test_pixscale_cd(): 496 mywcs = WCS(naxis=2) 497 mywcs.wcs.cd = [[-0.1, 0], [0, 0.2]] 498 mywcs.wcs.ctype = ['RA---TAN', 'DEC--TAN'] 499 assert_almost_equal(proj_plane_pixel_scales(mywcs), (0.1, 0.2)) 500 501 502@pytest.mark.parametrize('angle', 503 (30, 45, 60, 75)) 504def test_pixscale_cd_rotated(angle): 505 mywcs = WCS(naxis=2) 506 rho = np.radians(angle) 507 scale = 0.1 508 mywcs.wcs.cd = [[scale * np.cos(rho), -scale * np.sin(rho)], 509 [scale * np.sin(rho), scale * np.cos(rho)]] 510 mywcs.wcs.ctype = ['RA---TAN', 'DEC--TAN'] 511 assert_almost_equal(proj_plane_pixel_scales(mywcs), (0.1, 0.1)) 512 513 514@pytest.mark.parametrize('angle', 515 (30, 45, 60, 75)) 516def test_pixscale_pc_rotated(angle): 517 mywcs = WCS(naxis=2) 518 rho = np.radians(angle) 519 scale = 0.1 520 mywcs.wcs.cdelt = [-scale, scale] 521 mywcs.wcs.pc = [[np.cos(rho), -np.sin(rho)], 522 [np.sin(rho), np.cos(rho)]] 523 mywcs.wcs.ctype = ['RA---TAN', 'DEC--TAN'] 524 assert_almost_equal(proj_plane_pixel_scales(mywcs), (0.1, 0.1)) 525 526 527@pytest.mark.parametrize(('cdelt', 'pc', 'pccd'), 528 (([0.1, 0.2], np.eye(2), np.diag([0.1, 0.2])), 529 ([0.1, 0.2, 0.3], np.eye(3), np.diag([0.1, 0.2, 0.3])), 530 ([1, 1, 1], np.diag([0.1, 0.2, 0.3]), np.diag([0.1, 0.2, 0.3])))) 531def test_pixel_scale_matrix(cdelt, pc, pccd): 532 533 mywcs = WCS(naxis=(len(cdelt))) 534 mywcs.wcs.cdelt = cdelt 535 mywcs.wcs.pc = pc 536 537 assert_almost_equal(mywcs.pixel_scale_matrix, pccd) 538 539 540@pytest.mark.parametrize(('ctype', 'cel'), 541 ((['RA---TAN', 'DEC--TAN'], True), 542 (['RA---TAN', 'DEC--TAN', 'FREQ'], False), 543 (['RA---TAN', 'FREQ'], False),)) 544def test_is_celestial(ctype, cel): 545 mywcs = WCS(naxis=len(ctype)) 546 mywcs.wcs.ctype = ctype 547 548 assert mywcs.is_celestial == cel 549 550 551@pytest.mark.parametrize(('ctype', 'cel'), 552 ((['RA---TAN', 'DEC--TAN'], True), 553 (['RA---TAN', 'DEC--TAN', 'FREQ'], True), 554 (['RA---TAN', 'FREQ'], False),)) 555def test_has_celestial(ctype, cel): 556 mywcs = WCS(naxis=len(ctype)) 557 mywcs.wcs.ctype = ctype 558 559 assert mywcs.has_celestial == cel 560 561 562def test_has_celestial_correlated(): 563 # Regression test for astropy/astropy#8416 - has_celestial failed when 564 # celestial axes were correlated with other axes. 565 mywcs = WCS(naxis=3) 566 mywcs.wcs.ctype = 'RA---TAN', 'DEC--TAN', 'FREQ' 567 mywcs.wcs.cd = np.ones((3, 3)) 568 mywcs.wcs.set() 569 assert mywcs.has_celestial 570 571 572@pytest.mark.parametrize(('cdelt', 'pc', 'cd'), 573 ((np.array([0.1, 0.2]), np.eye(2), np.eye(2)), 574 (np.array([1, 1]), np.diag([0.1, 0.2]), np.eye(2)), 575 (np.array([0.1, 0.2]), np.eye(2), None), 576 (np.array([0.1, 0.2]), None, np.eye(2)), 577 )) 578def test_noncelestial_scale(cdelt, pc, cd): 579 580 mywcs = WCS(naxis=2) 581 if cd is not None: 582 mywcs.wcs.cd = cd 583 if pc is not None: 584 mywcs.wcs.pc = pc 585 586 # TODO: Some inputs emit RuntimeWarning from here onwards. 587 # Fix the test data. See @nden's comment in PR 9010. 588 with pytest.warns(None) as warning_lines: 589 mywcs.wcs.cdelt = cdelt 590 for w in warning_lines: 591 assert issubclass(w.category, RuntimeWarning) 592 assert 'cdelt will be ignored since cd is present' in str(w.message) 593 594 mywcs.wcs.ctype = ['RA---TAN', 'FREQ'] 595 596 ps = non_celestial_pixel_scales(mywcs) 597 598 assert_almost_equal(ps.to_value(u.deg), np.array([0.1, 0.2])) 599 600 601@pytest.mark.parametrize('mode', ['all', 'wcs']) 602def test_skycoord_to_pixel(mode): 603 604 # Import astropy.coordinates here to avoid circular imports 605 from astropy.coordinates import SkyCoord 606 607 header = get_pkg_data_contents('data/maps/1904-66_TAN.hdr', encoding='binary') 608 wcs = WCS(header) 609 610 ref = SkyCoord(0.1 * u.deg, -89. * u.deg, frame='icrs') 611 612 xp, yp = skycoord_to_pixel(ref, wcs, mode=mode) 613 614 # WCS is in FK5 so we need to transform back to ICRS 615 new = pixel_to_skycoord(xp, yp, wcs, mode=mode).transform_to('icrs') 616 617 assert_allclose(new.ra.degree, ref.ra.degree) 618 assert_allclose(new.dec.degree, ref.dec.degree) 619 620 # Make sure you can specify a different class using ``cls`` keyword 621 class SkyCoord2(SkyCoord): 622 pass 623 624 new2 = pixel_to_skycoord(xp, yp, wcs, mode=mode, 625 cls=SkyCoord2).transform_to('icrs') 626 627 assert new2.__class__ is SkyCoord2 628 assert_allclose(new2.ra.degree, ref.ra.degree) 629 assert_allclose(new2.dec.degree, ref.dec.degree) 630 631 632def test_skycoord_to_pixel_swapped(): 633 634 # Regression test for a bug that caused skycoord_to_pixel and 635 # pixel_to_skycoord to not work correctly if the axes were swapped in the 636 # WCS. 637 638 # Import astropy.coordinates here to avoid circular imports 639 from astropy.coordinates import SkyCoord 640 641 header = get_pkg_data_contents('data/maps/1904-66_TAN.hdr', encoding='binary') 642 wcs = WCS(header) 643 644 wcs_swapped = wcs.sub([WCSSUB_LATITUDE, WCSSUB_LONGITUDE]) 645 646 ref = SkyCoord(0.1 * u.deg, -89. * u.deg, frame='icrs') 647 648 xp1, yp1 = skycoord_to_pixel(ref, wcs) 649 xp2, yp2 = skycoord_to_pixel(ref, wcs_swapped) 650 651 assert_allclose(xp1, xp2) 652 assert_allclose(yp1, yp2) 653 654 # WCS is in FK5 so we need to transform back to ICRS 655 new1 = pixel_to_skycoord(xp1, yp1, wcs).transform_to('icrs') 656 new2 = pixel_to_skycoord(xp1, yp1, wcs_swapped).transform_to('icrs') 657 658 assert_allclose(new1.ra.degree, new2.ra.degree) 659 assert_allclose(new1.dec.degree, new2.dec.degree) 660 661 662def test_is_proj_plane_distorted(): 663 # non-orthogonal CD: 664 wcs = WCS(naxis=2) 665 wcs.wcs.cd = [[-0.1, 0], [0, 0.2]] 666 wcs.wcs.ctype = ['RA---TAN', 'DEC--TAN'] 667 assert(is_proj_plane_distorted(wcs)) 668 669 # almost orthogonal CD: 670 wcs.wcs.cd = [[0.1 + 2.0e-7, 1.7e-7], [1.2e-7, 0.1 - 1.3e-7]] 671 assert(not is_proj_plane_distorted(wcs)) 672 673 # real case: 674 header = get_pkg_data_filename('data/sip.fits') 675 with pytest.warns(FITSFixedWarning): 676 wcs = WCS(header) 677 assert(is_proj_plane_distorted(wcs)) 678 679 680@pytest.mark.parametrize('mode', ['all', 'wcs']) 681def test_skycoord_to_pixel_distortions(mode): 682 683 # Import astropy.coordinates here to avoid circular imports 684 from astropy.coordinates import SkyCoord 685 686 header = get_pkg_data_filename('data/sip.fits') 687 with pytest.warns(FITSFixedWarning): 688 wcs = WCS(header) 689 690 ref = SkyCoord(202.50 * u.deg, 47.19 * u.deg, frame='icrs') 691 692 xp, yp = skycoord_to_pixel(ref, wcs, mode=mode) 693 694 # WCS is in FK5 so we need to transform back to ICRS 695 new = pixel_to_skycoord(xp, yp, wcs, mode=mode).transform_to('icrs') 696 697 assert_allclose(new.ra.degree, ref.ra.degree) 698 assert_allclose(new.dec.degree, ref.dec.degree) 699 700 701@pytest.fixture 702def spatial_wcs_2d_small_angle(): 703 """ 704 This WCS has an almost linear correlation between the pixel and world axes 705 close to the reference pixel. 706 """ 707 wcs = WCS(naxis=2) 708 wcs.wcs.ctype = ['HPLN-TAN', 'HPLT-TAN'] 709 wcs.wcs.crpix = [3.0] * 2 710 wcs.wcs.cdelt = [0.002] * 2 711 wcs.wcs.crval = [0] * 2 712 wcs.wcs.set() 713 return wcs 714 715 716def test_local_pixel_derivatives(spatial_wcs_2d_small_angle): 717 not_diag = np.logical_not(np.diag([1, 1])) 718 # At (or close to) the reference pixel this should equal the cdelt 719 derivs = local_partial_pixel_derivatives(spatial_wcs_2d_small_angle, 3, 3) 720 np.testing.assert_allclose(np.diag(derivs), spatial_wcs_2d_small_angle.wcs.cdelt) 721 np.testing.assert_allclose(derivs[not_diag].flat, [0, 0], atol=1e-10) 722 723 # Far away from the reference pixel this should not equal the cdelt 724 derivs = local_partial_pixel_derivatives(spatial_wcs_2d_small_angle, 3e4, 3e4) 725 assert not np.allclose(np.diag(derivs), spatial_wcs_2d_small_angle.wcs.cdelt) 726 727 # At (or close to) the reference pixel this should equal the cdelt 728 derivs = local_partial_pixel_derivatives( 729 spatial_wcs_2d_small_angle, 3, 3, normalize_by_world=True) 730 np.testing.assert_allclose(np.diag(derivs), [1, 1]) 731 np.testing.assert_allclose(derivs[not_diag].flat, [0, 0], atol=1e-8) 732 733 734def test_pixel_to_world_correlation_matrix_celestial(): 735 736 wcs = WCS(naxis=2) 737 wcs.wcs.ctype = 'RA---TAN', 'DEC--TAN' 738 wcs.wcs.set() 739 740 assert_equal(wcs.axis_correlation_matrix, [[1, 1], [1, 1]]) 741 matrix, classes = _pixel_to_world_correlation_matrix(wcs) 742 assert_equal(matrix, [[1, 1]]) 743 assert classes == [SkyCoord] 744 745 746def test_pixel_to_world_correlation_matrix_spectral_cube_uncorrelated(): 747 748 wcs = WCS(naxis=3) 749 wcs.wcs.ctype = 'RA---TAN', 'FREQ', 'DEC--TAN' 750 wcs.wcs.set() 751 752 assert_equal(wcs.axis_correlation_matrix, [[1, 0, 1], [0, 1, 0], [1, 0, 1]]) 753 matrix, classes = _pixel_to_world_correlation_matrix(wcs) 754 assert_equal(matrix, [[1, 0, 1], [0, 1, 0]]) 755 assert classes == [SkyCoord, Quantity] 756 757 758def test_pixel_to_world_correlation_matrix_spectral_cube_correlated(): 759 760 wcs = WCS(naxis=3) 761 wcs.wcs.ctype = 'RA---TAN', 'FREQ', 'DEC--TAN' 762 wcs.wcs.cd = np.ones((3, 3)) 763 wcs.wcs.set() 764 765 assert_equal(wcs.axis_correlation_matrix, [[1, 1, 1], [1, 1, 1], [1, 1, 1]]) 766 matrix, classes = _pixel_to_world_correlation_matrix(wcs) 767 assert_equal(matrix, [[1, 1, 1], [1, 1, 1]]) 768 assert classes == [SkyCoord, Quantity] 769 770 771def test_pixel_to_pixel_correlation_matrix_celestial(): 772 773 wcs_in = WCS(naxis=2) 774 wcs_in.wcs.ctype = 'RA---TAN', 'DEC--TAN' 775 wcs_in.wcs.set() 776 777 wcs_out = WCS(naxis=2) 778 wcs_out.wcs.ctype = 'DEC--TAN', 'RA---TAN' 779 wcs_out.wcs.set() 780 781 matrix = _pixel_to_pixel_correlation_matrix(wcs_in, wcs_out) 782 assert_equal(matrix, [[1, 1], [1, 1]]) 783 784 785def test_pixel_to_pixel_correlation_matrix_spectral_cube_uncorrelated(): 786 787 wcs_in = WCS(naxis=3) 788 wcs_in.wcs.ctype = 'RA---TAN', 'DEC--TAN', 'FREQ' 789 wcs_in.wcs.set() 790 791 wcs_out = WCS(naxis=3) 792 wcs_out.wcs.ctype = 'DEC--TAN', 'FREQ', 'RA---TAN' 793 wcs_out.wcs.set() 794 795 matrix = _pixel_to_pixel_correlation_matrix(wcs_in, wcs_out) 796 assert_equal(matrix, [[1, 1, 0], [0, 0, 1], [1, 1, 0]]) 797 798 799def test_pixel_to_pixel_correlation_matrix_spectral_cube_correlated(): 800 801 # NOTE: only make one of the WCSes have correlated axes to really test this 802 803 wcs_in = WCS(naxis=3) 804 wcs_in.wcs.ctype = 'RA---TAN', 'DEC--TAN', 'FREQ' 805 wcs_in.wcs.set() 806 807 wcs_out = WCS(naxis=3) 808 wcs_out.wcs.ctype = 'DEC--TAN', 'FREQ', 'RA---TAN' 809 wcs_out.wcs.cd = np.ones((3, 3)) 810 wcs_out.wcs.set() 811 812 matrix = _pixel_to_pixel_correlation_matrix(wcs_in, wcs_out) 813 assert_equal(matrix, [[1, 1, 1], [1, 1, 1], [1, 1, 1]]) 814 815 816def test_pixel_to_pixel_correlation_matrix_mismatch(): 817 818 wcs_in = WCS(naxis=2) 819 wcs_in.wcs.ctype = 'RA---TAN', 'DEC--TAN' 820 wcs_in.wcs.set() 821 822 wcs_out = WCS(naxis=3) 823 wcs_out.wcs.ctype = 'DEC--TAN', 'FREQ', 'RA---TAN' 824 wcs_out.wcs.set() 825 826 with pytest.raises(ValueError) as exc: 827 _pixel_to_pixel_correlation_matrix(wcs_in, wcs_out) 828 assert exc.value.args[0] == "The two WCS return a different number of world coordinates" 829 830 wcs3 = WCS(naxis=2) 831 wcs3.wcs.ctype = 'FREQ', 'PIXEL' 832 wcs3.wcs.set() 833 834 with pytest.raises(ValueError) as exc: 835 _pixel_to_pixel_correlation_matrix(wcs_out, wcs3) 836 assert exc.value.args[0] == "The world coordinate types of the two WCS do not match" 837 838 wcs4 = WCS(naxis=4) 839 wcs4.wcs.ctype = 'RA---TAN', 'DEC--TAN', 'Q1', 'Q2' 840 wcs4.wcs.cunit = ['deg', 'deg', 'm/s', 'm/s'] 841 wcs4.wcs.set() 842 843 wcs5 = WCS(naxis=4) 844 wcs5.wcs.ctype = 'Q1', 'RA---TAN', 'DEC--TAN', 'Q2' 845 wcs5.wcs.cunit = ['m/s', 'deg', 'deg', 'm/s'] 846 wcs5.wcs.set() 847 848 with pytest.raises(ValueError, match="World coordinate order doesn't match " 849 "and automatic matching is ambiguous"): 850 _pixel_to_pixel_correlation_matrix(wcs4, wcs5) 851 852 853def test_pixel_to_pixel_correlation_matrix_nonsquare(): 854 855 # Here we set up an input WCS that maps 3 pixel coordinates to 4 world 856 # coordinates - the idea is to make sure that things work fine in cases 857 # where the number of input and output pixel coordinates do not match. 858 859 class FakeWCS(object): 860 pass 861 862 wcs_in = FakeWCS() 863 wcs_in.low_level_wcs = wcs_in 864 wcs_in.pixel_n_dim = 3 865 wcs_in.world_n_dim = 4 866 wcs_in.axis_correlation_matrix = [[True, True, False], 867 [True, True, False], 868 [True, True, False], 869 [False, False, True]] 870 wcs_in.world_axis_object_components = [('spat', 'ra', 'ra.degree'), 871 ('spat', 'dec', 'dec.degree'), 872 ('spec', 0, 'value'), 873 ('time', 0, 'utc.value')] 874 wcs_in.world_axis_object_classes = {'spat': ('astropy.coordinates.SkyCoord', (), 875 {'frame': 'icrs'}), 876 'spec': ('astropy.units.Wavelength', (None,), {}), 877 'time': ('astropy.time.Time', (None,), 878 {'format': 'mjd', 'scale': 'utc'})} 879 880 wcs_out = FakeWCS() 881 wcs_out.low_level_wcs = wcs_out 882 wcs_out.pixel_n_dim = 4 883 wcs_out.world_n_dim = 4 884 wcs_out.axis_correlation_matrix = [[True, False, False, False], 885 [False, True, True, False], 886 [False, True, True, False], 887 [False, False, False, True]] 888 wcs_out.world_axis_object_components = [('spec', 0, 'value'), 889 ('spat', 'ra', 'ra.degree'), 890 ('spat', 'dec', 'dec.degree'), 891 ('time', 0, 'utc.value')] 892 wcs_out.world_axis_object_classes = wcs_in.world_axis_object_classes 893 894 matrix = _pixel_to_pixel_correlation_matrix(wcs_in, wcs_out) 895 896 matrix = matrix.astype(int) 897 898 # The shape should be (n_pixel_out, n_pixel_in) 899 assert matrix.shape == (4, 3) 900 901 expected = np.array([[1, 1, 0], [1, 1, 0], [1, 1, 0], [0, 0, 1]]) 902 assert_equal(matrix, expected) 903 904 905def test_split_matrix(): 906 907 assert _split_matrix(np.array([[1]])) == [([0], [0])] 908 909 assert _split_matrix(np.array([[1, 1], 910 [1, 1]])) == [([0, 1], [0, 1])] 911 912 assert _split_matrix(np.array([[1, 1, 0], 913 [1, 1, 0], 914 [0, 0, 1]])) == [([0, 1], [0, 1]), ([2], [2])] 915 916 assert _split_matrix(np.array([[0, 1, 0], 917 [1, 0, 0], 918 [0, 0, 1]])) == [([0], [1]), ([1], [0]), ([2], [2])] 919 920 assert _split_matrix(np.array([[0, 1, 1], 921 [1, 0, 0], 922 [1, 0, 1]])) == [([0, 1, 2], [0, 1, 2])] 923 924 925def test_pixel_to_pixel(): 926 927 wcs_in = WCS(naxis=3) 928 wcs_in.wcs.ctype = 'DEC--TAN', 'FREQ', 'RA---TAN' 929 wcs_in.wcs.set() 930 931 wcs_out = WCS(naxis=3) 932 wcs_out.wcs.ctype = 'GLON-CAR', 'GLAT-CAR', 'FREQ' 933 wcs_out.wcs.set() 934 935 # First try with scalars 936 with pytest.warns(AstropyUserWarning, match='No observer defined on WCS'): 937 x, y, z = pixel_to_pixel(wcs_in, wcs_out, 1, 2, 3) 938 assert x.shape == () 939 assert y.shape == () 940 assert z.shape == () 941 942 # Now try with broadcasted arrays 943 x = np.linspace(10, 20, 10) 944 y = np.linspace(10, 20, 20) 945 z = np.linspace(10, 20, 30) 946 Z1, Y1, X1 = np.meshgrid(z, y, x, indexing='ij', copy=False) 947 with pytest.warns(AstropyUserWarning, match='No observer defined on WCS'): 948 X2, Y2, Z2 = pixel_to_pixel(wcs_in, wcs_out, X1, Y1, Z1) 949 950 # The final arrays should have the correct shape 951 assert X2.shape == (30, 20, 10) 952 assert Y2.shape == (30, 20, 10) 953 assert Z2.shape == (30, 20, 10) 954 955 # But behind the scenes should also be broadcasted 956 assert unbroadcast(X2).shape == (30, 1, 10) 957 assert unbroadcast(Y2).shape == (30, 1, 10) 958 assert unbroadcast(Z2).shape == (20, 1) 959 960 # We can put the values back through the function to ensure round-tripping 961 with pytest.warns(AstropyUserWarning, match='No observer defined on WCS'): 962 X3, Y3, Z3 = pixel_to_pixel(wcs_out, wcs_in, X2, Y2, Z2) 963 964 # The final arrays should have the correct shape 965 assert X2.shape == (30, 20, 10) 966 assert Y2.shape == (30, 20, 10) 967 assert Z2.shape == (30, 20, 10) 968 969 # But behind the scenes should also be broadcasted 970 assert unbroadcast(X3).shape == (30, 1, 10) 971 assert unbroadcast(Y3).shape == (20, 1) 972 assert unbroadcast(Z3).shape == (30, 1, 10) 973 974 # And these arrays should match the input 975 assert_allclose(X1, X3) 976 assert_allclose(Y1, Y3) 977 assert_allclose(Z1, Z3) 978 979 980def test_pixel_to_pixel_correlated(): 981 982 wcs_in = WCS(naxis=2) 983 wcs_in.wcs.ctype = 'DEC--TAN', 'RA---TAN' 984 wcs_in.wcs.set() 985 986 wcs_out = WCS(naxis=2) 987 wcs_out.wcs.ctype = 'GLON-CAR', 'GLAT-CAR' 988 wcs_out.wcs.set() 989 990 # First try with scalars 991 x, y = pixel_to_pixel(wcs_in, wcs_out, 1, 2) 992 assert x.shape == () 993 assert y.shape == () 994 995 # Now try with broadcasted arrays 996 x = np.linspace(10, 20, 10) 997 y = np.linspace(10, 20, 20) 998 Y1, X1 = np.meshgrid(y, x, indexing='ij', copy=False) 999 Y2, X2 = pixel_to_pixel(wcs_in, wcs_out, X1, Y1) 1000 1001 # The final arrays should have the correct shape 1002 assert X2.shape == (20, 10) 1003 assert Y2.shape == (20, 10) 1004 1005 # and there are no efficiency gains here since the celestial axes are correlated 1006 assert unbroadcast(X2).shape == (20, 10) 1007 1008 1009def test_pixel_to_pixel_1d(): 1010 1011 # Simple test to make sure that when WCS only returns one world coordinate 1012 # this still works correctly (since this requires special treatment behind 1013 # the scenes). 1014 1015 wcs_in = WCS(naxis=1) 1016 wcs_in.wcs.ctype = 'COORD1', 1017 wcs_in.wcs.cunit = 'nm', 1018 wcs_in.wcs.set() 1019 1020 wcs_out = WCS(naxis=1) 1021 wcs_out.wcs.ctype = 'COORD2', 1022 wcs_out.wcs.cunit = 'cm', 1023 wcs_out.wcs.set() 1024 1025 # First try with a scalar 1026 x = pixel_to_pixel(wcs_in, wcs_out, 1) 1027 assert x.shape == () 1028 1029 # Next with a regular array 1030 x = np.linspace(10, 20, 10) 1031 x = pixel_to_pixel(wcs_in, wcs_out, x) 1032 assert x.shape == (10,) 1033 1034 # And now try with a broadcasted array 1035 x = np.broadcast_to(np.linspace(10, 20, 10), (4, 10)) 1036 x = pixel_to_pixel(wcs_in, wcs_out, x) 1037 assert x.shape == (4, 10) 1038 1039 # The broadcasting of the input should be retained 1040 assert unbroadcast(x).shape == (10,) 1041 1042 1043header_str_linear = """ 1044XTENSION= 'IMAGE ' / Image extension 1045BITPIX = -32 / array data type 1046NAXIS = 2 / number of array dimensions 1047NAXIS1 = 50 1048NAXIS2 = 50 1049PCOUNT = 0 / number of parameters 1050GCOUNT = 1 / number of groups 1051RADESYS = 'ICRS ' 1052EQUINOX = 2000.0 1053WCSAXES = 2 1054CTYPE1 = 'RA---TAN' 1055CTYPE2 = 'DEC--TAN' 1056CRVAL1 = 250.3497414839765 1057CRVAL2 = 2.280925599609063 1058CRPIX1 = 1045.0 1059CRPIX2 = 1001.0 1060CD1_1 = -0.005564478186178 1061CD1_2 = -0.001042099258152 1062CD2_1 = 0.00118144146585 1063CD2_2 = -0.005590816683583 1064""" 1065header_str_sip = """ 1066XTENSION= 'IMAGE ' / Image extension 1067BITPIX = -32 / array data type 1068NAXIS = 2 / number of array dimensions 1069NAXIS1 = 50 1070NAXIS2 = 50 1071PCOUNT = 0 / number of parameters 1072GCOUNT = 1 / number of groups 1073RADESYS = 'ICRS ' 1074EQUINOX = 2000.0 1075WCSAXES = 2 1076CTYPE1 = 'RA---TAN-SIP' 1077CTYPE2 = 'DEC--TAN-SIP' 1078CRVAL1 = 250.3497414839765 1079CRVAL2 = 2.280925599609063 1080CRPIX1 = 1045.0 1081CRPIX2 = 1001.0 1082CD1_1 = -0.005564478186178 1083CD1_2 = -0.001042099258152 1084CD2_1 = 0.00118144146585 1085CD2_2 = -0.005590816683583 1086A_ORDER = 2 1087B_ORDER = 2 1088A_2_0 = 2.02451189234E-05 1089A_0_2 = 3.317603337918E-06 1090A_1_1 = 1.73456334971071E-05 1091B_2_0 = 3.331330003472E-06 1092B_0_2 = 2.04247482482589E-05 1093B_1_1 = 1.71476710804143E-05 1094AP_ORDER= 2 1095BP_ORDER= 2 1096AP_1_0 = 0.000904700296389636 1097AP_0_1 = 0.000627660715584716 1098AP_2_0 = -2.023482905861E-05 1099AP_0_2 = -3.332285841011E-06 1100AP_1_1 = -1.731636633824E-05 1101BP_1_0 = 0.000627960882053211 1102BP_0_1 = 0.000911222886084808 1103BP_2_0 = -3.343918167224E-06 1104BP_0_2 = -2.041598249021E-05 1105BP_1_1 = -1.711876336719E-05 1106A_DMAX = 44.72893589844534 1107B_DMAX = 44.62692873032506 1108""" 1109header_str_prob = """ 1110NAXIS = 2 / number of array dimensions 1111WCSAXES = 2 / Number of coordinate axes 1112CRPIX1 = 1024.5 / Pixel coordinate of reference point 1113CRPIX2 = 1024.5 / Pixel coordinate of reference point 1114CD1_1 = -1.7445934400771E-05 / Coordinate transformation matrix element 1115CD1_2 = -4.9826985362578E-08 / Coordinate transformation matrix element 1116CD2_1 = -5.0068838822312E-08 / Coordinate transformation matrix element 1117CD2_2 = 1.7530614610951E-05 / Coordinate transformation matrix element 1118CTYPE1 = 'RA---TAN' / Right ascension, gnomonic projection 1119CTYPE2 = 'DEC--TAN' / Declination, gnomonic projection 1120CRVAL1 = 5.8689341666667 / [deg] Coordinate value at reference point 1121CRVAL2 = -71.995508583333 / [deg] Coordinate value at reference point 1122""" 1123 1124 1125@pytest.mark.skipif('not HAS_SCIPY') 1126@pytest.mark.parametrize( 1127 'header_str,crval,sip_degree,user_proj_point,exp_max_dist,exp_std_dist', 1128 [ 1129 # simple testset no distortions 1130 (header_str_linear, 250.3497414839765, None, False, 7e-5*u.deg, 2.5e-5*u.deg), 1131 # simple testset with distortions 1132 (header_str_sip, 250.3497414839765, 2, False, 7e-6*u.deg, 2.5e-6*u.deg), 1133 # testset with problematic WCS header that failed before 1134 (header_str_prob, 5.8689341666667, None, False, 7e-6*u.deg, 2.5e-6*u.deg), 1135 # simple testset no distortions, user defined center 1136 (header_str_linear, 250.3497414839765, None, True, 7e-5*u.deg, 2.5e-5*u.deg), 1137 # 360->0 degree crossover, simple testset no distortions 1138 (header_str_linear, 352.3497414839765, None, False, 7e-5*u.deg, 2.5e-5*u.deg), 1139 # 360->0 degree crossover, simple testset with distortions 1140 (header_str_sip, 352.3497414839765, 2, False, 7e-6*u.deg, 2.5e-6*u.deg), 1141 # 360->0 degree crossover, testset with problematic WCS header that failed before 1142 (header_str_prob, 352.3497414839765, None, False, 7e-6*u.deg, 2.5e-6*u.deg), 1143 # 360->0 degree crossover, simple testset no distortions, user defined center 1144 (header_str_linear, 352.3497414839765, None, True, 7e-5*u.deg, 2.5e-5*u.deg), 1145 ]) 1146def test_fit_wcs_from_points(header_str, crval, sip_degree, user_proj_point, 1147 exp_max_dist, exp_std_dist): 1148 header = fits.Header.fromstring(header_str, sep='\n') 1149 header["CRVAL1"] = crval 1150 1151 true_wcs = WCS(header, relax=True) 1152 1153 # Getting the pixel coordinates 1154 x, y = np.meshgrid(list(range(10)), list(range(10))) 1155 x = x.flatten() 1156 y = y.flatten() 1157 1158 # Calculating the true sky positions 1159 world_pix = true_wcs.pixel_to_world(x, y) 1160 1161 # which projection point to use 1162 if user_proj_point: 1163 proj_point = world_pix[0] 1164 projlon = proj_point.data.lon.deg 1165 projlat = proj_point.data.lat.deg 1166 else: 1167 proj_point = 'center' 1168 1169 # Fitting the wcs 1170 fit_wcs = fit_wcs_from_points((x, y), world_pix, 1171 proj_point=proj_point, 1172 sip_degree=sip_degree) 1173 1174 # Validate that the true sky coordinates 1175 # match sky coordinates calculated from the wcs fit 1176 world_pix_new = fit_wcs.pixel_to_world(x, y) 1177 1178 dists = world_pix.separation(world_pix_new) 1179 1180 assert dists.max() < exp_max_dist 1181 assert np.std(dists) < exp_std_dist 1182 1183 if user_proj_point: 1184 assert (fit_wcs.wcs.crval == [projlon, projlat]).all() 1185 1186 1187@pytest.mark.skipif('not HAS_SCIPY') 1188def test_fit_wcs_from_points_CRPIX_bounds(): 1189 # Test CRPIX bounds requirement 1190 wcs_str = """ 1191WCSAXES = 2 / Number of coordinate axes 1192CRPIX1 = 1045.0 / Pixel coordinate of reference point 1193CRPIX2 = 1001.0 / Pixel coordinate of reference point 1194PC1_1 = 0.00056205870415378 / Coordinate transformation matrix element 1195PC1_2 = -0.00569181083243 / Coordinate transformation matrix element 1196PC2_1 = 0.0056776810932466 / Coordinate transformation matrix element 1197PC2_2 = 0.0004208048403273 / Coordinate transformation matrix element 1198CDELT1 = 1.0 / [deg] Coordinate increment at reference point 1199CDELT2 = 1.0 / [deg] Coordinate increment at reference point 1200CUNIT1 = 'deg' / Units of coordinate increment and value 1201CUNIT2 = 'deg' / Units of coordinate increment and value 1202CTYPE1 = 'RA---TAN' / Right ascension, gnomonic projection 1203CTYPE2 = 'DEC--TAN' / Declination, gnomonic projection 1204CRVAL1 = 104.57797893504 / [deg] Coordinate value at reference point 1205CRVAL2 = -74.195502593322 / [deg] Coordinate value at reference point 1206LONPOLE = 180.0 / [deg] Native longitude of celestial pole 1207LATPOLE = -74.195502593322 / [deg] Native latitude of celestial pole 1208TIMESYS = 'TDB' / Time scale 1209TIMEUNIT= 'd' / Time units 1210DATEREF = '1858-11-17' / ISO-8601 fiducial time 1211MJDREFI = 0.0 / [d] MJD of fiducial time, integer part 1212MJDREFF = 0.0 / [d] MJD of fiducial time, fractional part 1213DATE-OBS= '2019-03-27T03:30:13.832Z' / ISO-8601 time of observation 1214MJD-OBS = 58569.145993426 / [d] MJD of observation 1215MJD-OBS = 58569.145993426 / [d] MJD at start of observation 1216TSTART = 1569.6467941661 / [d] Time elapsed since fiducial time at start 1217DATE-END= '2019-03-27T04:00:13.831Z' / ISO-8601 time at end of observation 1218MJD-END = 58569.166826748 / [d] MJD at end of observation 1219TSTOP = 1569.6676274905 / [d] Time elapsed since fiducial time at end 1220TELAPSE = 0.02083332443 / [d] Elapsed time (start to stop) 1221TIMEDEL = 0.020833333333333 / [d] Time resolution 1222TIMEPIXR= 0.5 / Reference position of timestamp in binned data 1223RADESYS = 'ICRS' / Equatorial coordinate system 1224""" 1225 wcs_header = fits.Header.fromstring(wcs_str, sep='\n') 1226 ffi_wcs = WCS(wcs_header) 1227 1228 yi, xi = (1000, 1000) 1229 y, x = (10, 200) 1230 1231 center_coord = SkyCoord(ffi_wcs.all_pix2world([[xi+x//2, yi+y//2]], 0), unit='deg')[0] 1232 ypix, xpix = [arr.flatten() for arr in np.mgrid[xi : xi + x, yi : yi + y]] 1233 world_pix = SkyCoord(*ffi_wcs.all_pix2world(xpix, ypix, 0), unit='deg') 1234 1235 fit_wcs = fit_wcs_from_points((ypix, xpix), world_pix, proj_point='center') 1236 1237 assert (fit_wcs.wcs.crpix.astype(int) == [1100, 1005]).all() 1238 assert fit_wcs.pixel_shape == (1199, 1009) 1239 1240 1241@pytest.mark.skipif('not HAS_SCIPY') 1242def test_issue10991(): 1243 # test issue #10991 (it just needs to run and set the user defined crval) 1244 xy = np.array([[1766.88276168, 662.96432257, 171.50212526, 120.70924648], 1245 [1706.69832901, 1788.85480559, 1216.98949653, 1307.41843381]]) 1246 world_coords = SkyCoord([(66.3542367 , 22.20000162), (67.15416174, 19.18042906), 1247 (65.73375432, 17.54251555), (66.02400512, 17.44413253)], 1248 frame="icrs", unit="deg") 1249 proj_point = SkyCoord(64.67514918, 19.63389538, 1250 frame="icrs", unit="deg") 1251 1252 fit_wcs = fit_wcs_from_points( 1253 xy=xy, 1254 world_coords=world_coords, 1255 proj_point=proj_point, 1256 projection='TAN' 1257 ) 1258 projlon = proj_point.data.lon.deg 1259 projlat = proj_point.data.lat.deg 1260 assert (fit_wcs.wcs.crval == [projlon, projlat]).all() 1261 1262 1263@pytest.mark.remote_data 1264@pytest.mark.parametrize('x_in,y_in', [[0, 0], [np.arange(5), np.arange(5)]]) 1265def test_pixel_to_world_itrs(x_in, y_in): 1266 """Regression test for https://github.com/astropy/astropy/pull/9609""" 1267 with pytest.warns(None) as w: 1268 wcs = WCS({'NAXIS': 2, 1269 'CTYPE1': 'TLON-CAR', 1270 'CTYPE2': 'TLAT-CAR', 1271 'RADESYS': 'ITRS ', 1272 'DATE-OBS': '2017-08-17T12:41:04.444'}) 1273 1274 if Version(_wcs.__version__) >= Version('7.4'): 1275 assert len(w) == 1 1276 msg = str(w[0].message) 1277 assert "'datfix' made the change 'Set MJD-OBS to 57982.528524 from DATE-OBS'." in msg 1278 else: 1279 assert len(w) == 0 1280 1281 # This shouldn't raise an exception. 1282 coord = wcs.pixel_to_world(x_in, y_in) 1283 1284 # Check round trip transformation. 1285 x, y = wcs.world_to_pixel(coord) 1286 1287 np.testing.assert_almost_equal(x, x_in) 1288 np.testing.assert_almost_equal(y, y_in) 1289 1290 1291@pytest.fixture 1292def dkist_location(): 1293 return EarthLocation(*(-5466045.25695494, -2404388.73741278, 2242133.88769004) * u.m) 1294 1295 1296def test_obsgeo_cartesian(dkist_location): 1297 obstime = Time("2021-05-21T03:00:00") 1298 wcs = WCS(naxis=2) 1299 wcs.wcs.obsgeo = list(dkist_location.to_value(u.m).tolist()) + [0, 0, 0] 1300 wcs.wcs.dateobs = obstime.isot 1301 1302 frame = obsgeo_to_frame(wcs.wcs.obsgeo, obstime) 1303 1304 assert isinstance(frame, ITRS) 1305 assert frame.x == dkist_location.x 1306 assert frame.y == dkist_location.y 1307 assert frame.z == dkist_location.z 1308 1309 1310def test_obsgeo_spherical(dkist_location): 1311 obstime = Time("2021-05-21T03:00:00") 1312 dkist_location = dkist_location.get_itrs(obstime) 1313 loc_sph = dkist_location.spherical 1314 1315 wcs = WCS(naxis=2) 1316 wcs.wcs.obsgeo = [0, 0, 0] + [loc_sph.lon.value, loc_sph.lat.value, loc_sph.distance.value] 1317 wcs.wcs.dateobs = obstime.isot 1318 1319 frame = obsgeo_to_frame(wcs.wcs.obsgeo, obstime) 1320 1321 assert isinstance(frame, ITRS) 1322 assert u.allclose(frame.x, dkist_location.x) 1323 assert u.allclose(frame.y, dkist_location.y) 1324 assert u.allclose(frame.z, dkist_location.z) 1325 1326 1327def test_obsgeo_infinite(dkist_location): 1328 obstime = Time("2021-05-21T03:00:00") 1329 dkist_location = dkist_location.get_itrs(obstime) 1330 loc_sph = dkist_location.spherical 1331 1332 wcs = WCS(naxis=2) 1333 wcs.wcs.obsgeo = [1, 1, np.nan] + [loc_sph.lon.value, loc_sph.lat.value, loc_sph.distance.value] 1334 wcs.wcs.dateobs = obstime.isot 1335 wcs.wcs.set() 1336 1337 frame = obsgeo_to_frame(wcs.wcs.obsgeo, obstime) 1338 1339 assert isinstance(frame, ITRS) 1340 assert u.allclose(frame.x, dkist_location.x) 1341 assert u.allclose(frame.y, dkist_location.y) 1342 assert u.allclose(frame.z, dkist_location.z) 1343 1344 1345@pytest.mark.parametrize("obsgeo", ([np.nan] * 6, None, [0] * 6, [54] * 5)) 1346def test_obsgeo_invalid(obsgeo): 1347 1348 with pytest.raises(ValueError): 1349 obsgeo_to_frame(obsgeo, None) 1350