DM-39168: Add utility to compute apparent motion #65
Conversation
cmsaunders
force-pushed
the
tickets/DM-39168
branch
3 times, most recently
from
bcce819 to
af0fde4
Compare
There was a problem hiding this comment.
Just to check: why can't you use astropy to do the PM+parallax correction directly? I see you're doing just that in the tests; would it make sense to just make this function a wrapper around astropy.SkyCoord conversions?
I'll trust you on the parallax math itself: I don't have the reference to hand.
tests/test_gbdesAstrometricFit.py
Outdated
| astropy_dDec = newCoords.dec.degree - refNewCoords.dec.degree | ||
|
|
||
| np.testing.assert_allclose(astropy_dRA, ra_correction.value, rtol=1e-1) | ||
| np.testing.assert_allclose(astropy_dDec, dec_correction.value, rtol=1e-1) |
There was a problem hiding this comment.
1e-1 degrees accuracy doesn't seem very good?
Separately, the mathematical definition for rtol/atol in np.allclose always trips me up. I'd suggest using our own assertFloatsAlmostEqual, as the rtol/atol distinction is clearer (it's one or the other, not a combination of the two).
There was a problem hiding this comment.
Yeah, this is why I asked you to review the ticket! They're also both different from the value from gbdes, which is done via a change of coordinate systems. If I recall correctly, if you start digging into the astropy calculation, it takes you into a maze of ERFA code, and I never figured out exactly how the parallax part was being calculated. Maybe @timj has a better idea for how to approach this?
For the allclose vs assertFloatsAlmostEqual, I just forgot that assertFloatsAlmostEqual allows arrays, unlike assertAlmostEqual, so I have switched to that.
There was a problem hiding this comment.
The upshot is that the actual differences are all about -0.012 mas, which may come down to the choice of ecliptic angle between the function and astropy, and I have decided this is acceptable in this context.
There was a problem hiding this comment.
I am still worrying that your Times are all in UTC and would like confirmation that they really should be in UTC and not TAI.
There was a problem hiding this comment.
I'm still replying to comments and haven't pushed any changes yet. The times come from the VisitTable and should be TAI. This is endemic in gbdesAstrometricFit.
There was a problem hiding this comment.
Ok. So every astropy.time.Time you create will need scale="tai" added to it.
There was a problem hiding this comment.
Thanks, changing to assertFloatsAlmostEqual should make the exact comparison more obvious. And Tim's right about requiring scale in the astropy time objects.
There was a problem hiding this comment.
Right, I already made the scale change.
|
It's not possible to use the |
tests/test_gbdesAstrometricFit.py
Outdated
| dec=self.Decs * u.degree, | ||
| pm_ra_cosdec=pmRA_cosDec * u.mas / u.yr, | ||
| pm_dec=self.pmDecs * u.mas / u.yr, | ||
| obstime=Time(refEpoch, format="mjd"), |
There was a problem hiding this comment.
General point -- be careful with timescales. If you know you are using tai everywhere you might want to specify that. If you know it's UTC somewhere, then say so. All your Time objects are defaulting to UTC and that might not be what you mean given that all our times in FITS metadata and elsewhere are TAI.
cmsaunders
force-pushed
the
tickets/DM-39168
branch
2 times, most recently
from
e7c1634 to
bf23a25
Compare
tests/test_gbdesAstrometricFit.py
Outdated
| # Set parallaxes to zero to test proper motion part by itself. | ||
| parallaxes = np.zeros(len(self.ras)) * u.mas | ||
|
|
||
| data = Table( |
There was a problem hiding this comment.
I would recommend importing this with the full namespace as well, as Table is another generic term like Time, otherwise it's not obvious that it came from astropy.
| astropyDRA = newCoords.ra.degree - coords.ra.degree | ||
| astropyDDec = newCoords.dec.degree - coords.dec.degree | ||
|
|
||
| self.assertFloatsAlmostEqual(astropyDRA, raCorrection.value, rtol=1e-6) |
There was a problem hiding this comment.
Maybe make this comparison directly in milliarcsecond, as I can't translate 1e-6 into mas in my head to know whether this is a "good" result or not. Something like this should make it more obvious how much the difference is, though what exactly you need for atol I don't know, but I think that's what you really want here.
| self.assertFloatsAlmostEqual(astropyDRA, raCorrection.value, rtol=1e-6) | |
| self.assertFloatsAlmostEqual(astropyDRA.to_value(u.mas), raCorrection.to_value(u.mas), atol=1) |
There was a problem hiding this comment.
This is the relative difference though, so 1e-6 is unitless. I added a comment that the proper motions are on the scale of 0-3 mas/yr and the differences are on the scale of 1e-8 mas.
tests/test_gbdesAstrometricFit.py
Outdated
| astropyDRA = newCoords.ra.degree - refNewCoords.ra.degree | ||
| astropyDDec = newCoords.dec.degree - refNewCoords.dec.degree | ||
|
|
||
| self.assertFloatsAlmostEqual(astropyDRA, raCorrection.value, rtol=1e-1) |
There was a problem hiding this comment.
As I noted above, you should do the comparison in mas and probably with atol, so we can see how accurate it is in "real" units. It's not obvious to me what a 10% accuracy means here, but it doesn't feel like it's very good?
There was a problem hiding this comment.
Well, the actual accuracy is about 1.1% at worst, but I rounded up rather than have the value be too tuned. I guess I could do 2e-2. For the atol, I'm not sure that it's preferable here, since you don't know what the values of astropyDRA and raCorrection are here (they are ~1-2 mas). I'm going to just add a comment on the scale of the errors here, which correspond to 0.04 mas at most.
There was a problem hiding this comment.
Ah, thank you for that clarification. That's not as bad as I'd thought.
cmsaunders
force-pushed
the
tickets/DM-39168
branch
from
bf23a25 to
26b0488
Compare
No description provided.