get_constellation() should use exact HMS, DMS boundaries, not 4-decimal approximations

[nealmcb]

Description

IAU's official constellation boundaries are defined using HMS / DMS coordinates in
E. Delporte | Délimitation scientifique des constellations, cartes, Cambridge University Press, London 1930.

But the data that drives the astropy get_constellation() routine is derived by Roman et al. (1987) from Delporte's work, using a lossy conversion to decimal hourangles and degrees with just 4 decimal places.

Roman's data is archived at http://cdsarc.u-strasbg.fr/viz-bin/Cat?VI/42
astropy's copy is in constellation_data_roman87.dat file.

For example, the 4th line is:

astropy/astropy/coordinates/data/constellation_data_roman87.dat

Line 4 in 51c619e

21.0000 23.0000 86.1667 UMi

It represents an edge between Ursa Minor (UMi there) and Cepheus.

It matches data on page 14 of Delporte's book, as reproduced at http://www.atlascoelestis.com/del%20inice%20cost.htm, where
the corresponding data shows up as the 2nd line after the Ursa Minor sequence starts:

Paralléle de 86° 10' de 21 h. o m. à 23 h. o m.

We thus see that the original parallel at dec 86° 10' was rounded off by Roman to 86.1667.

(Note that edges are listed twice in Delporte's book, once for each constellation, so it also shows up on page 39 for UMi.)

A better digital representation of Delporte's IAU constellation standard is provided in multiple formats at Pierre Barbier's Constellation Boundaries

For example, line 195 of his edges_18.txt file shows the same edge represented exactly in the native sexagesimal format:

318:319 P+ 21:00:00 +86:10:00 23:00:00 +86:10:00 CEP UMI

Based on the roundoff, astropy thus would find the wrong constellation for coordinates like 22h 86.16668°

While trying to demonstrate this, I ran across what seems to be another issue that I understand less well. get_constellation() calls PrecessedGeocentric(equinox='B1875') which throws the coordinates off a bit as I document below, interpreting them as geocentric with a distance inferred somehow I guess.

astropy/astropy/coordinates/funcs.py

Line 235 in b1addf5

constel_coord = coord.transform_to(PrecessedGeocentric(equinox='B1875'))

While this might make sense for a SkyCoord pointing to a comet or other solar system object, it makes no sense to me when trying to figure out which constellation a distant variable star is in. Perhaps there is a more appropriate way to call the functions, but the methods I'm using below are similar to those in test_constellations() in test_funcs.py, so I'm still confused. Note that those tests don't test the edge case (pun intended) that I'm discussing there. I have some additional tests which will, which I'll add to that file in a PR.

To be precise:

Actual behavior

from astropy.coordinates import SkyCoord, FK5, get_constellation, PrecessedGeocentric
from astropy import units as u

print('Coordinate\tTransformed\tConstellation\tHand annotation')

for dec in [86.16, 86.1666, 86.16668, 86.1668]:
    c = FK5(22*u.hour, dec*u.deg, equinox='B1875')
    constel = get_constellation(c)
    pg = c.transform_to(PrecessedGeocentric(equinox='B1875'))
    print(f'{dec:<10.8}\t{pg.dec:.8}\t{constel}')

Annotated output:

Coordinate      Transformed     Constellation   Hand annotation
86.16           86.164512 deg   Cepheus     # Right
86.1666         86.171112 deg   Ursa Minor  # Wrong
86.16668        86.171192 deg   Ursa Minor  # Right for wrong reason
86.1668         86.171312 deg   Ursa Minor  # Right

Expected behavior

I expected the 2nd row (86.1666) to properly be located in Cepheus.
I expected the transformed column (showing the values that are actually looked up in the table, currently generated when PrecessedGeocentric has been applied to the coordinate) above to match the Coordinate column.
After the issue with PrecessedGeocentric is resolved, I expect the 3rd row to continue to show up in Ursa Minor.

Steps to Reproduce

See code above.

Notes

I also note that at http://pbarbier.com/terms.html Barbier says "All the original content on his site is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License, so we should be able to re-use it.

System Details

Linux-4.15.0-72-generic-x86_64-with-Ubuntu-18.04-bionic
Python 3.6.9 (default, Nov 7 2019, 10:44:02)
[GCC 8.3.0]
Numpy 1.18.1
astropy 4.0

[nealmcb]

To fix the second issue I ran into, would it work to just replace

astropy/astropy/coordinates/funcs.py

Line 235 in b1addf5

constel_coord = coord.transform_to(PrecessedGeocentric(equinox='B1875'))

with

constel_coord = coord.transform_to(FK5(equinox='B1875'))

Offhand that doesn't seem to cause the geocentrism(?) problems I noted above.

And with that I get the answers I was expecting from the badly transformed Roman data.

Coordinate      Constellation
86.16           Cepheus
86.1666         Cepheus
86.16668        Cepheus
86.1668         Ursa Minor

I.e. the answer for 86.16668 is wrong, for the reasons I started off clarifying.
The right answers are Cepheus for the first two, and Ursa Minor for the second two.

[eteq]

I think your approach is right, @nealmcb ! I think I'm the one who did the PrecessedGeocentric, and that was based on an assumption I think you've revealed here is incorrect about how to interpret the Roman paper - i.e., I guessed that this was meant to be interpreted in the geocentric way because it's based on where given coordinates "were" in B1875. But I hadn't appreciated that Delporte's work was the "true" answer, so I think your proposed change is correct based on the above.

[eteq]

This was incorrectly auto-closed in #9856 because #9856 documents the poroblem, but doesn't solve it. So next step is to actually implement the plan @nealmcb suggested above to make the #9856 no longer xfail

[nealmcb]

Thanks for accepting the failing test!

I've explored the bugs and possible fixes a bit more in a Jupyter notebook, which I've put online for clarity:
Explore astropy get_constellation bugs and Table ASCII CDS format

This is preliminary, but it seems that the bug with the improper use of PrecessedGeocentric affects perhaps 0.03% of a random set of queries, many-fold more than the 0.002% of queries affected by the original bug I identified with rounding off the boundary values.

Also, it seems that (as explained at Constellation a celestial body is in) the binary search Python code in SkyField (API Reference) is much faster than the sequential approach in astropy (based on Roman). See also Fast Determination of Constellation Membership - NASA/ADS. So we might want to just build off that code (MIT license, it seems).

I'm curious to query some published catalogs with constellation information to see if any of them reflect errors derived from astropy. Perhaps constellation lookups are rare enough that this problem hasn't made much difference, but it seems possible.

[nealmcb]

Here is a quick sense of the performance difference between the astropy get_constellation() code and SkyField's binary search-based constellation_at() code. The latter is over 100x faster in this quick comparison via the IPython %%timeit magic function!

ra, dec = randomRaDec()
p = FK5(ra*u.hour, dec*u.deg, equinox='B1875')
%%timeit
get_constellation(p, short_name=True)
4.89 ms ± 76.4 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

%%timeit
ra, dec = randomRaDec()
constellation_at(position_of_radec(ra, dec, epoch=B1875))
37.3 µs ± 853 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)