Fix for erfa 1.7.2 by mhvk · Pull Request #11485 · astropy/astropy

mhvk · 2021-04-06T01:29:05Z

This fixes most of the problems with ERFA 1.7.2, following the suggestion of @mkbrewer to adjust erfa_astrom.apio to follow the update in ERFA itself. With that, we're internally consistent.

I'm a bit puzzled, though, by the remaining failure, which is a comparison against pyephem. Note that the actual failure is against an earlier astropy result, but why is it off by so much? Also comparing directly to the pyephem result, we're also off by a large amount (490"), much larger than used to be the case.

~~Now I noticed one strange thing with that particular test, which is that the height is 30000 m. However, a quick check with 30 m did not help. (link to test below)~~

EDIT: as noted by @StuartLittlefair below, the really strange number is the input altitude of -60 degrees! With that more reasonable, things pass.

Also a question: should I make apio ERFA-version dependent, or can we simply insist on a recent ERFA version (I think definitely OK for 4.3, but what about LTS?).

See

astropy/astropy/coordinates/tests/accuracy/test_altaz_icrs.py

Lines 94 to 127 in fff5ccf

    
           def test_against_pyephem(): 
        
               """Check that Astropy gives consistent results with one PyEphem example. 
        
               PyEphem: https://rhodesmill.org/pyephem/ 
        
               See example input and output here: 
        
               https://gist.github.com/zonca/1672906 
        
               https://github.com/phn/pytpm/issues/2#issuecomment-3698679 
        
               """ 
        
               obstime = Time('2011-09-18 08:50:00') 
        
               location = EarthLocation(lon=Angle('-109d24m53.1s'), 
        
                                        lat=Angle('33d41m46.0s'), 
        
                                        height=30000. * u.m) 
        
               # We are using the default pressure and temperature in PyEphem 
        
               # relative_humidity = ? 
        
               # obswl = ? 
        
               altaz_frame = AltAz(obstime=obstime, location=location, 
        
                                   temperature=15 * u.deg_C, pressure=1.010 * u.bar) 
        
               altaz = SkyCoord('6.8927d -60.7665d', frame=altaz_frame) 
        
               radec_actual = altaz.transform_to('icrs') 
        
               radec_expected = SkyCoord('196.497518d -4.569323d', frame='icrs')  # EPHEM 
        
               # radec_expected = SkyCoord('196.496220d -4.569390d', frame='icrs')  # HORIZON 
        
               distance = radec_actual.separation(radec_expected).to('arcsec') 
        
               # TODO: why is this difference so large? 
        
               # It currently is: 31.45187984720655 arcsec 
        
               assert distance < 1e3 * u.arcsec 
        
               # Add assert on current Astropy result so that we notice if something changes 
        
               radec_expected = SkyCoord('196.495372d -4.560694d', frame='icrs') 
        
               distance = radec_actual.separation(radec_expected).to('arcsec') 
        
               # Current value: 0.0031402822944751997 arcsec 
        
               assert distance < 1 * u.arcsec

github-actions · 2021-04-06T01:29:59Z

👋 Thank you for your draft pull request! Do you know that you can use [ci skip] or [skip ci] in your commit messages to skip running continuous integration tests until you are ready?

StuartLittlefair · 2021-04-06T10:03:14Z

I'm very confused by this myself, not least how this test passed in previous releases. I'm on astropy 4.2.1 and it is 489.6" away from the pypehem value and 521.28" away from the "current" astropy result.

BTW if you take the "expected" RA/Dec from PyEphem and transform it into the altaz frame with astropy you get

<SkyCoord (AltAz: obstime=2011-09-18 08:50:00.000, location=(-1774037.41919622, -5033515.95039486, 3535095.77328894) m, pressure=1010.0 hPa, temperature=15.0 deg_C, relative_humidity=0.0, obswl=1.0 micron): (az, alt) in deg
    (6.89010736, -60.74182357)>

If you do the same thing with SkyField you get az=6.890133 and alt=-60.7666 degrees. So whilst there isn't perfect agreement, both astropy and Skyfield agree that the azimuth required to get the "expected" RA/Dec from PyEphem would need to be 6.8901, not 6.8927.

We seem to have a bigger problem with altitude. Skyfield gets this about right, but we are quite far off. Things aren't improved in astropy master c.f 4.2.1 (master uses the direct ICRS->AltAz transform).

However, I still haven't understood why

pllim · 2021-04-06T13:52:34Z

@adrn has a weekly coordinates tag-up now. Do you think it would be helpful to bring this issue up in the next tag-up?

adrn · 2021-04-06T13:57:03Z

Next week we're doing the triaging meeting (right?), so I added this to the agenda for Apr 19!

StuartLittlefair · 2021-04-06T14:03:19Z

It is worth discussing, but I think the issue is to do with an abuse of the atmospheric refraction code. It isn't really designed to work at altitudes of -60 degrees! If you turn off atmospheric refraction the disagreement drops to 4.5 arc seconds!

StuartLittlefair · 2021-04-06T14:08:40Z

By the way I'm ok with the changes in this PR, and happy to approve, but I think we should clear up why this test is failing now and how it passed in prior releases before merging.

pllim · 2021-04-06T14:09:24Z

Next week we're doing the triaging meeting (right?)

Ah, yes. I was gonna get in touch with you about it today. Let's take this part offline.

mhvk · 2021-04-06T19:15:23Z

Hmm, I don't think we have to expect to reproduce altitudes well below the horizon! Surely that is a mistake! (But let me ping @cdeil who put in this test). I think I'll just run pyephem for a more normal altitude and put that in...

mhvk · 2021-04-06T19:37:11Z

OK, I made a commit where I changed the pyephem input parameters to something a bit more sensible.

mhvk · 2021-04-06T19:50:43Z


-
-@pytest.mark.parametrize('icoo', icrs_coords)
-def test_icrs_consistency(icoo):


I have no idea why the test as written iterated over all the SkyCoord rather than just do them in one go. Changes testing time from ~10 s to <1 s.

StuartLittlefair · 2021-04-07T10:53:50Z

So I've approved this PR because the code is absolutely fine. With regards to whether we can require pyerfa>=1.7.2 I personally do not see any problems this could cause with LTS. However, maybe we should ask someone who actually has to package these things, like @olebole ?

mhvk · 2021-04-07T14:14:15Z

Note that I can quite easily change it such that the code becomes erfa-version dependent (like we do for numpy).

olebole · 2021-04-08T13:47:02Z

@StuartLittlefair Debian is currently in Hard Freeze for the new release 11 "Bullseye"; therefore I will not upload anything unless it fixes a release critical bug. And for our new release cycle, requiring >= 1.7.2 is absolutely fine -- 1.7.2 is what we will have in Bullseye, so it would even not prevent backporting.

mhvk · 2021-04-11T20:56:24Z

OK, I think for 4.3 we should just require erfa >= 1.7.2. I had worried about 4.0.x, but since erfa_astrom didn't exist in 4.0, that would need a different solution anyway. EDIT: if there is even a problem, since we didn't replace it with our faster version.

Since this was approved before, I'll merge once tests pass again - that's helpful too for getting pyerfa 1.7.3 out - this time I'll test it works with astropy before releasing!

mhvk added coordinates erfa labels Apr 6, 2021

mhvk added this to the v4.3 milestone Apr 6, 2021

mhvk requested a review from StuartLittlefair April 6, 2021 01:29

github-actions Bot added the installation label Apr 6, 2021

mhvk mentioned this pull request Apr 6, 2021

test_icrs_consistency and test_against_pyephem failing with pyerfa 1.7.2 #11305

Closed

1 task done

mhvk commented Apr 6, 2021

View reviewed changes

StuartLittlefair approved these changes Apr 7, 2021

View reviewed changes

mhvk marked this pull request as ready for review April 11, 2021 20:48

mhvk added 5 commits April 11, 2021 16:54

Remove limitation on pyerfa

d7b4ed0

Adjust apio following erfa update.

fdd69b2

Fix bad precision test against pyephem

c827dd2

Do all ICRS to observed tests in one go.

4b8a49c

Add changelog fragment notifying the minimum requirement

7d9da0f

mhvk force-pushed the fix-for-erfa-1.7.2 branch from 3a6839e to 7d9da0f Compare April 11, 2021 20:54

mhvk merged commit 522d5cd into astropy:main Apr 11, 2021

mhvk deleted the fix-for-erfa-1.7.2 branch April 11, 2021 23:16

	def test_against_pyephem():
	"""Check that Astropy gives consistent results with one PyEphem example.

	PyEphem: https://rhodesmill.org/pyephem/

	See example input and output here:
	https://gist.github.com/zonca/1672906
	https://github.com/phn/pytpm/issues/2#issuecomment-3698679
	"""
	obstime = Time('2011-09-18 08:50:00')
	location = EarthLocation(lon=Angle('-109d24m53.1s'),
	lat=Angle('33d41m46.0s'),
	height=30000. * u.m)
	# We are using the default pressure and temperature in PyEphem
	# relative_humidity = ?
	# obswl = ?
	altaz_frame = AltAz(obstime=obstime, location=location,
	temperature=15 * u.deg_C, pressure=1.010 * u.bar)

	altaz = SkyCoord('6.8927d -60.7665d', frame=altaz_frame)
	radec_actual = altaz.transform_to('icrs')

	radec_expected = SkyCoord('196.497518d -4.569323d', frame='icrs') # EPHEM
	# radec_expected = SkyCoord('196.496220d -4.569390d', frame='icrs') # HORIZON
	distance = radec_actual.separation(radec_expected).to('arcsec')
	# TODO: why is this difference so large?
	# It currently is: 31.45187984720655 arcsec
	assert distance < 1e3 * u.arcsec

	# Add assert on current Astropy result so that we notice if something changes
	radec_expected = SkyCoord('196.495372d -4.560694d', frame='icrs')
	distance = radec_actual.separation(radec_expected).to('arcsec')
	# Current value: 0.0031402822944751997 arcsec
	assert distance < 1 * u.arcsec



		@pytest.mark.parametrize('icoo', icrs_coords)
		def test_icrs_consistency(icoo):

Conversation

mhvk commented Apr 6, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

github-actions Bot commented Apr 6, 2021

Uh oh!

StuartLittlefair commented Apr 6, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

pllim commented Apr 6, 2021

Uh oh!

adrn commented Apr 6, 2021

Uh oh!

StuartLittlefair commented Apr 6, 2021

Uh oh!

StuartLittlefair commented Apr 6, 2021

Uh oh!

pllim commented Apr 6, 2021

Uh oh!

mhvk commented Apr 6, 2021

Uh oh!

mhvk commented Apr 6, 2021

Uh oh!

mhvk Apr 6, 2021

Choose a reason for hiding this comment

Uh oh!

StuartLittlefair commented Apr 7, 2021

Uh oh!

mhvk commented Apr 7, 2021

Uh oh!

olebole commented Apr 8, 2021

Uh oh!

mhvk commented Apr 11, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

5 participants

mhvk commented Apr 6, 2021 •

edited

Loading

StuartLittlefair commented Apr 6, 2021 •

edited

Loading

mhvk commented Apr 11, 2021 •

edited

Loading