1# -*- coding: utf-8 -*- 2# Licensed under a 3-clause BSD style license - see LICENSE.rst 3 4""" 5This is the APE5 coordinates API document re-written to work as a series of test 6functions. 7 8Note that new tests for coordinates functionality should generally *not* be 9added to this file - instead, add them to other appropriate test modules in 10this package, like ``test_sky_coord.py``, ``test_frames.py``, or 11``test_representation.py``. This file is instead meant mainly to keep track of 12deviations from the original APE5 plan. 13""" 14 15import pytest 16import numpy as np 17from numpy import testing as npt 18 19from astropy.tests.helper import assert_quantity_allclose as assert_allclose 20from astropy import units as u 21from astropy import time 22from astropy import coordinates as coords 23from astropy.units import allclose 24from astropy.utils.compat.optional_deps import HAS_SCIPY # noqa 25 26 27def test_representations_api(): 28 from astropy.coordinates.representation import SphericalRepresentation, \ 29 UnitSphericalRepresentation, PhysicsSphericalRepresentation, \ 30 CartesianRepresentation 31 from astropy.coordinates import Angle, Longitude, Latitude, Distance 32 33 # <-----------------Classes for representation of coordinate data--------------> 34 # These classes inherit from a common base class and internally contain Quantity 35 # objects, which are arrays (although they may act as scalars, like numpy's 36 # length-0 "arrays") 37 38 # They can be initialized with a variety of ways that make intuitive sense. 39 # Distance is optional. 40 UnitSphericalRepresentation(lon=8*u.hour, lat=5*u.deg) 41 UnitSphericalRepresentation(lon=8*u.hourangle, lat=5*u.deg) 42 SphericalRepresentation(lon=8*u.hourangle, lat=5*u.deg, distance=10*u.kpc) 43 44 # In the initial implementation, the lat/lon/distance arguments to the 45 # initializer must be in order. A *possible* future change will be to allow 46 # smarter guessing of the order. E.g. `Latitude` and `Longitude` objects can be 47 # given in any order. 48 UnitSphericalRepresentation(Longitude(8, u.hour), Latitude(5, u.deg)) 49 SphericalRepresentation(Longitude(8, u.hour), Latitude(5, u.deg), Distance(10, u.kpc)) 50 51 # Arrays of any of the inputs are fine 52 UnitSphericalRepresentation(lon=[8, 9]*u.hourangle, lat=[5, 6]*u.deg) 53 54 # Default is to copy arrays, but optionally, it can be a reference 55 UnitSphericalRepresentation(lon=[8, 9]*u.hourangle, lat=[5, 6]*u.deg, copy=False) 56 57 # strings are parsed by `Latitude` and `Longitude` constructors, so no need to 58 # implement parsing in the Representation classes 59 UnitSphericalRepresentation(lon=Angle('2h6m3.3s'), lat=Angle('0.1rad')) 60 61 # Or, you can give `Quantity`s with keywords, and they will be internally 62 # converted to Angle/Distance 63 c1 = SphericalRepresentation(lon=8*u.hourangle, lat=5*u.deg, distance=10*u.kpc) 64 65 # Can also give another representation object with the `reprobj` keyword. 66 c2 = SphericalRepresentation.from_representation(c1) 67 68 # distance, lat, and lon typically will just match in shape 69 SphericalRepresentation(lon=[8, 9]*u.hourangle, lat=[5, 6]*u.deg, distance=[10, 11]*u.kpc) 70 # if the inputs are not the same, if possible they will be broadcast following 71 # numpy's standard broadcasting rules. 72 c2 = SphericalRepresentation(lon=[8, 9]*u.hourangle, lat=[5, 6]*u.deg, distance=10*u.kpc) 73 assert len(c2.distance) == 2 74 # when they can't be broadcast, it is a ValueError (same as Numpy) 75 with pytest.raises(ValueError): 76 c2 = UnitSphericalRepresentation(lon=[8, 9, 10]*u.hourangle, lat=[5, 6]*u.deg) 77 78 # It's also possible to pass in scalar quantity lists with mixed units. These 79 # are converted to array quantities following the same rule as `Quantity`: all 80 # elements are converted to match the first element's units. 81 c2 = UnitSphericalRepresentation(lon=Angle([8*u.hourangle, 135*u.deg]), 82 lat=Angle([5*u.deg, (6*np.pi/180)*u.rad])) 83 assert c2.lat.unit == u.deg and c2.lon.unit == u.hourangle 84 npt.assert_almost_equal(c2.lon[1].value, 9) 85 86 # The Quantity initializer itself can also be used to force the unit even if the 87 # first element doesn't have the right unit 88 lon = u.Quantity([120*u.deg, 135*u.deg], u.hourangle) 89 lat = u.Quantity([(5*np.pi/180)*u.rad, 0.4*u.hourangle], u.deg) 90 c2 = UnitSphericalRepresentation(lon, lat) 91 92 # regardless of how input, the `lat` and `lon` come out as angle/distance 93 assert isinstance(c1.lat, Angle) 94 assert isinstance(c1.lat, Latitude) # `Latitude` is an `~astropy.coordinates.Angle` subclass 95 assert isinstance(c1.distance, Distance) 96 97 # but they are read-only, as representations are immutable once created 98 with pytest.raises(AttributeError): 99 c1.lat = Latitude(5, u.deg) 100 # Note that it is still possible to modify the array in-place, but this is not 101 # sanctioned by the API, as this would prevent things like caching. 102 c2.lat[:] = [0] * u.deg # possible, but NOT SUPPORTED 103 104 # To address the fact that there are various other conventions for how spherical 105 # coordinates are defined, other conventions can be included as new classes. 106 # Later there may be other conventions that we implement - for now just the 107 # physics convention, as it is one of the most common cases. 108 _ = PhysicsSphericalRepresentation(phi=120*u.deg, theta=85*u.deg, r=3*u.kpc) 109 110 # first dimension must be length-3 if a lone `Quantity` is passed in. 111 c1 = CartesianRepresentation(np.random.randn(3, 100) * u.kpc) 112 assert c1.xyz.shape[0] == 3 113 assert c1.xyz.unit == u.kpc 114 assert c1.x.shape[0] == 100 115 assert c1.y.shape[0] == 100 116 assert c1.z.shape[0] == 100 117 # can also give each as separate keywords 118 CartesianRepresentation(x=np.random.randn(100)*u.kpc, 119 y=np.random.randn(100)*u.kpc, 120 z=np.random.randn(100)*u.kpc) 121 # if the units don't match but are all distances, they will automatically be 122 # converted to match `x` 123 xarr, yarr, zarr = np.random.randn(3, 100) 124 c1 = CartesianRepresentation(x=xarr*u.kpc, y=yarr*u.kpc, z=zarr*u.kpc) 125 c2 = CartesianRepresentation(x=xarr*u.kpc, y=yarr*u.kpc, z=zarr*u.pc) 126 assert c1.xyz.unit == c2.xyz.unit == u.kpc 127 assert_allclose((c1.z / 1000) - c2.z, 0*u.kpc, atol=1e-10*u.kpc) 128 129 # representations convert into other representations via `represent_as` 130 srep = SphericalRepresentation(lon=90*u.deg, lat=0*u.deg, distance=1*u.pc) 131 crep = srep.represent_as(CartesianRepresentation) 132 assert_allclose(crep.x, 0*u.pc, atol=1e-10*u.pc) 133 assert_allclose(crep.y, 1*u.pc, atol=1e-10*u.pc) 134 assert_allclose(crep.z, 0*u.pc, atol=1e-10*u.pc) 135 # The functions that actually do the conversion are defined via methods on the 136 # representation classes. This may later be expanded into a full registerable 137 # transform graph like the coordinate frames, but initially it will be a simpler 138 # method system 139 140 141def test_frame_api(): 142 from astropy.coordinates.representation import SphericalRepresentation, \ 143 UnitSphericalRepresentation 144 from astropy.coordinates.builtin_frames import ICRS, FK5 145 # <--------------------Reference Frame/"Low-level" classes---------------------> 146 # The low-level classes have a dual role: they act as specifiers of coordinate 147 # frames and they *may* also contain data as one of the representation objects, 148 # in which case they are the actual coordinate objects themselves. 149 150 # They can always accept a representation as a first argument 151 icrs = ICRS(UnitSphericalRepresentation(lon=8*u.hour, lat=5*u.deg)) 152 153 # which is stored as the `data` attribute 154 assert icrs.data.lat == 5*u.deg 155 assert icrs.data.lon == 8*u.hourangle 156 157 # Frames that require additional information like equinoxs or obstimes get them 158 # as keyword parameters to the frame constructor. Where sensible, defaults are 159 # used. E.g., FK5 is almost always J2000 equinox 160 fk5 = FK5(UnitSphericalRepresentation(lon=8*u.hour, lat=5*u.deg)) 161 J2000 = time.Time('J2000') 162 fk5_2000 = FK5(UnitSphericalRepresentation(lon=8*u.hour, lat=5*u.deg), equinox=J2000) 163 assert fk5.equinox == fk5_2000.equinox 164 165 # the information required to specify the frame is immutable 166 J2001 = time.Time('J2001') 167 with pytest.raises(AttributeError): 168 fk5.equinox = J2001 169 170 # Similar for the representation data. 171 with pytest.raises(AttributeError): 172 fk5.data = UnitSphericalRepresentation(lon=8*u.hour, lat=5*u.deg) 173 174 # There is also a class-level attribute that lists the attributes needed to 175 # identify the frame. These include attributes like `equinox` shown above. 176 assert all(nm in ('equinox', 'obstime') for nm in fk5.get_frame_attr_names()) 177 178 # the result of `get_frame_attr_names` is called for particularly in the 179 # high-level class (discussed below) to allow round-tripping between various 180 # frames. It is also part of the public API for other similar developer / 181 # advanced users' use. 182 183 # The actual position information is accessed via the representation objects 184 assert_allclose(icrs.represent_as(SphericalRepresentation).lat, 5*u.deg) 185 # shorthand for the above 186 assert_allclose(icrs.spherical.lat, 5*u.deg) 187 assert icrs.cartesian.z.value > 0 188 189 # Many frames have a "default" representation, the one in which they are 190 # conventionally described, often with a special name for some of the 191 # coordinates. E.g., most equatorial coordinate systems are spherical with RA and 192 # Dec. This works simply as a shorthand for the longer form above 193 194 assert_allclose(icrs.dec, 5*u.deg) 195 assert_allclose(fk5.ra, 8*u.hourangle) 196 197 assert icrs.representation_type == SphericalRepresentation 198 199 # low-level classes can also be initialized with names valid for that representation 200 # and frame: 201 icrs_2 = ICRS(ra=8*u.hour, dec=5*u.deg, distance=1*u.kpc) 202 assert_allclose(icrs.ra, icrs_2.ra) 203 204 # and these are taken as the default if keywords are not given: 205 # icrs_nokwarg = ICRS(8*u.hour, 5*u.deg, distance=1*u.kpc) 206 # assert icrs_nokwarg.ra == icrs_2.ra and icrs_nokwarg.dec == icrs_2.dec 207 208 # they also are capable of computing on-sky or 3d separations from each other, 209 # which will be a direct port of the existing methods: 210 coo1 = ICRS(ra=0*u.hour, dec=0*u.deg) 211 coo2 = ICRS(ra=0*u.hour, dec=1*u.deg) 212 # `separation` is the on-sky separation 213 assert_allclose(coo1.separation(coo2).degree, 1.0) 214 215 # while `separation_3d` includes the 3D distance information 216 coo3 = ICRS(ra=0*u.hour, dec=0*u.deg, distance=1*u.kpc) 217 coo4 = ICRS(ra=0*u.hour, dec=0*u.deg, distance=2*u.kpc) 218 assert coo3.separation_3d(coo4).kpc == 1.0 219 220 # The next example fails because `coo1` and `coo2` don't have distances 221 with pytest.raises(ValueError): 222 assert coo1.separation_3d(coo2).kpc == 1.0 223 224 # repr/str also shows info, with frame and data 225 # assert repr(fk5) == '' 226 227 228def test_transform_api(): 229 from astropy.coordinates.representation import UnitSphericalRepresentation 230 from astropy.coordinates.builtin_frames import ICRS, FK5 231 from astropy.coordinates.baseframe import frame_transform_graph, BaseCoordinateFrame 232 from astropy.coordinates.transformations import DynamicMatrixTransform 233 # <------------------------Transformations-------------------------------------> 234 # Transformation functionality is the key to the whole scheme: they transform 235 # low-level classes from one frame to another. 236 237 # (used below but defined above in the API) 238 fk5 = FK5(ra=8*u.hour, dec=5*u.deg) 239 240 # If no data (or `None`) is given, the class acts as a specifier of a frame, but 241 # without any stored data. 242 J2001 = time.Time('J2001') 243 fk5_J2001_frame = FK5(equinox=J2001) 244 245 # if they do not have data, the string instead is the frame specification 246 assert repr(fk5_J2001_frame) == "<FK5 Frame (equinox=J2001.000)>" 247 248 # Note that, although a frame object is immutable and can't have data added, it 249 # can be used to create a new object that does have data by giving the 250 # `realize_frame` method a representation: 251 srep = UnitSphericalRepresentation(lon=8*u.hour, lat=5*u.deg) 252 fk5_j2001_with_data = fk5_J2001_frame.realize_frame(srep) 253 assert fk5_j2001_with_data.data is not None 254 # Now `fk5_j2001_with_data` is in the same frame as `fk5_J2001_frame`, but it 255 # is an actual low-level coordinate, rather than a frame without data. 256 257 # These frames are primarily useful for specifying what a coordinate should be 258 # transformed *into*, as they are used by the `transform_to` method 259 # E.g., this snippet precesses the point to the new equinox 260 newfk5 = fk5.transform_to(fk5_J2001_frame) 261 assert newfk5.equinox == J2001 262 263 # transforming to a new frame necessarily loses framespec information if that 264 # information is not applicable to the new frame. This means transforms are not 265 # always round-trippable: 266 fk5_2 = FK5(ra=8*u.hour, dec=5*u.deg, equinox=J2001) 267 ic_trans = fk5_2.transform_to(ICRS()) 268 269 # `ic_trans` does not have an `equinox`, so now when we transform back to FK5, 270 # it's a *different* RA and Dec 271 fk5_trans = ic_trans.transform_to(FK5()) 272 assert not allclose(fk5_2.ra, fk5_trans.ra, rtol=0, atol=1e-10*u.deg) 273 274 # But if you explicitly give the right equinox, all is fine 275 fk5_trans_2 = fk5_2.transform_to(FK5(equinox=J2001)) 276 assert_allclose(fk5_2.ra, fk5_trans_2.ra, rtol=0, atol=1e-10*u.deg) 277 278 # Trying to transforming a frame with no data is of course an error: 279 with pytest.raises(ValueError): 280 FK5(equinox=J2001).transform_to(ICRS()) 281 282 # To actually define a new transformation, the same scheme as in the 283 # 0.2/0.3 coordinates framework can be re-used - a graph of transform functions 284 # connecting various coordinate classes together. The main changes are: 285 # 1) The transform functions now get the frame object they are transforming the 286 # current data into. 287 # 2) Frames with additional information need to have a way to transform between 288 # objects of the same class, but with different framespecinfo values 289 290 # An example transform function: 291 class SomeNewSystem(BaseCoordinateFrame): 292 pass 293 294 @frame_transform_graph.transform(DynamicMatrixTransform, SomeNewSystem, FK5) 295 def new_to_fk5(newobj, fk5frame): 296 _ = newobj.obstime 297 _ = fk5frame.equinox 298 # ... build a *cartesian* transform matrix using `eq` that transforms from 299 # the `newobj` frame as observed at `ot` to FK5 an equinox `eq` 300 matrix = np.eye(3) 301 return matrix 302 303 # Other options for transform functions include one that simply returns the new 304 # coordinate object, and one that returns a cartesian matrix but does *not* 305 # require `newobj` or `fk5frame` - this allows optimization of the transform. 306 307 308def test_highlevel_api(): 309 J2001 = time.Time('J2001') 310 311 # <--------------------------"High-level" class--------------------------------> 312 # The "high-level" class is intended to wrap the lower-level classes in such a 313 # way that they can be round-tripped, as well as providing a variety of 314 # convenience functionality. This document is not intended to show *all* of the 315 # possible high-level functionality, rather how the high-level classes are 316 # initialized and interact with the low-level classes 317 318 # this creates an object that contains an `ICRS` low-level class, initialized 319 # identically to the first ICRS example further up. 320 321 sc = coords.SkyCoord(coords.SphericalRepresentation(lon=8 * u.hour, 322 lat=5 * u.deg, distance=1 * u.kpc), frame='icrs') 323 324 # Other representations and `system` keywords delegate to the appropriate 325 # low-level class. The already-existing registry for user-defined coordinates 326 # will be used by `SkyCoordinate` to figure out what various the `system` 327 # keyword actually means. 328 329 sc = coords.SkyCoord(ra=8 * u.hour, dec=5 * u.deg, frame='icrs') 330 sc = coords.SkyCoord(l=120 * u.deg, b=5 * u.deg, frame='galactic') 331 332 # High-level classes can also be initialized directly from low-level objects 333 sc = coords.SkyCoord(coords.ICRS(ra=8 * u.hour, dec=5 * u.deg)) 334 335 # The next example raises an error because the high-level class must always 336 # have position data. 337 with pytest.raises(ValueError): 338 sc = coords.SkyCoord(coords.FK5(equinox=J2001)) # raises ValueError 339 340 # similarly, the low-level object can always be accessed 341 342 # this is how it's supposed to look, but sometimes the numbers get rounded in 343 # funny ways 344 # assert repr(sc.frame) == '<ICRS Coordinate: ra=120.0 deg, dec=5.0 deg>' 345 rscf = repr(sc.frame) 346 assert rscf.startswith('<ICRS Coordinate: (ra, dec) in deg') 347 348 # and the string representation will be inherited from the low-level class. 349 350 # same deal, should loook like this, but different archituectures/ python 351 # versions may round the numbers differently 352 # assert repr(sc) == '<SkyCoord (ICRS): ra=120.0 deg, dec=5.0 deg>' 353 rsc = repr(sc) 354 assert rsc.startswith('<SkyCoord (ICRS): (ra, dec) in deg') 355 356 # Supports a variety of possible complex string formats 357 sc = coords.SkyCoord('8h00m00s +5d00m00.0s', frame='icrs') 358 359 # In the next example, the unit is only needed b/c units are ambiguous. In 360 # general, we *never* accept ambiguity 361 sc = coords.SkyCoord('8:00:00 +5:00:00.0', unit=(u.hour, u.deg), frame='icrs') 362 363 # The next one would yield length-2 array coordinates, because of the comma 364 365 sc = coords.SkyCoord(['8h 5d', '2°2′3″ 0.3rad'], frame='icrs') 366 367 # It should also interpret common designation styles as a coordinate 368 # NOT YET 369 # sc = coords.SkyCoord('SDSS J123456.89-012345.6', frame='icrs') 370 371 # but it should also be possible to provide formats for outputting to strings, 372 # similar to `Time`. This can be added right away or at a later date. 373 374 # transformation is done the same as for low-level classes, which it delegates to 375 376 sc_fk5_j2001 = sc.transform_to(coords.FK5(equinox=J2001)) 377 assert sc_fk5_j2001.equinox == J2001 378 379 # The key difference is that the high-level class remembers frame information 380 # necessary for round-tripping, unlike the low-level classes: 381 sc1 = coords.SkyCoord(ra=8 * u.hour, dec=5 * u.deg, equinox=J2001, frame='fk5') 382 sc2 = sc1.transform_to('icrs') 383 384 # The next assertion succeeds, but it doesn't mean anything for ICRS, as ICRS 385 # isn't defined in terms of an equinox 386 assert sc2.equinox == J2001 387 388 # But it *is* necessary once we transform to FK5 389 sc3 = sc2.transform_to('fk5') 390 assert sc3.equinox == J2001 391 assert_allclose(sc1.ra, sc3.ra) 392 393 # `SkyCoord` will also include the attribute-style access that is in the 394 # v0.2/0.3 coordinate objects. This will *not* be in the low-level classes 395 sc = coords.SkyCoord(ra=8 * u.hour, dec=5 * u.deg, frame='icrs') 396 scgal = sc.galactic 397 assert str(scgal).startswith('<SkyCoord (Galactic): (l, b)') 398 399 # the existing `from_name` and `match_to_catalog_*` methods will be moved to the 400 # high-level class as convenience functionality. 401 402 # in remote-data test below! 403 # m31icrs = coords.SkyCoord.from_name('M31', frame='icrs') 404 # assert str(m31icrs) == '<SkyCoord (ICRS) RA=10.68471 deg, Dec=41.26875 deg>' 405 406 if HAS_SCIPY: 407 cat1 = coords.SkyCoord(ra=[1, 2]*u.hr, dec=[3, 4.01]*u.deg, 408 distance=[5, 6]*u.kpc, frame='icrs') 409 cat2 = coords.SkyCoord(ra=[1, 2, 2.01]*u.hr, dec=[3, 4, 5]*u.deg, 410 distance=[5, 200, 6]*u.kpc, frame='icrs') 411 idx1, sep2d1, dist3d1 = cat1.match_to_catalog_sky(cat2) 412 idx2, sep2d2, dist3d2 = cat1.match_to_catalog_3d(cat2) 413 414 assert np.any(idx1 != idx2) 415 416 # additional convenience functionality for the future should be added as methods 417 # on `SkyCoord`, *not* the low-level classes. 418 419 420@pytest.mark.remote_data 421def test_highlevel_api_remote(): 422 m31icrs = coords.SkyCoord.from_name('M31', frame='icrs') 423 424 m31str = str(m31icrs) 425 assert m31str.startswith('<SkyCoord (ICRS): (ra, dec) in deg\n (') 426 assert m31str.endswith(')>') 427 assert '10.68' in m31str 428 assert '41.26' in m31str 429 # The above is essentially a replacement of the below, but tweaked so that 430 # small/moderate changes in what `from_name` returns don't cause the tests 431 # to fail 432 # assert str(m31icrs) == '<SkyCoord (ICRS): (ra, dec) in deg\n (10.6847083, 41.26875)>' 433 434 m31fk4 = coords.SkyCoord.from_name('M31', frame='fk4') 435 436 assert not m31icrs.is_equivalent_frame(m31fk4) 437 assert np.abs(m31icrs.ra - m31fk4.ra) > .5*u.deg 438