Repository navigation
Improved light-deflection code can have issues for the center of the Sun #10878
Description
#10666 improved the handling of light deflection for coordinate transformation, but it apparently can have issues when transforming the center of the Sun. I've seen it happen once on the SunPy CI, where the RuntimeWarning is treated as an error. However, I haven't succeeded yet with reproducing it locally, presumably because it depends on a very precise selection of time (sounds like a job for hypothesis...). Here's an example traceback (failed build), where:
sun = SkyCoord(0*u.deg, 0*u.deg, 0*u.AU, frame='hcrs', obstime=time)
time.utc = "2020-10-15 17:21:50.208105"Traceback:
../../.tox/py38-devdeps/lib/python3.8/site-packages/sunpy/coordinates/sun.py:199: in apparent_longitude
coord = sun.transform_to(GeocentricMeanEcliptic(equinox=time))
../../.tox/py38-devdeps/lib/python3.8/site-packages/astropy/coordinates/sky_coordinate.py:619: in transform_to
new_coord = trans(self.frame, generic_frame)
../../.tox/py38-devdeps/lib/python3.8/site-packages/astropy/coordinates/transformations.py:1392: in __call__
curr_coord = t(curr_coord, curr_toframe)
../../.tox/py38-devdeps/lib/python3.8/site-packages/astropy/coordinates/transformations.py:995: in __call__
return supcall(fromcoord, toframe)
../../.tox/py38-devdeps/lib/python3.8/site-packages/astropy/coordinates/builtin_frames/icrs_cirs_transforms.py:135: in icrs_to_gcrs
gcrs_ra, gcrs_dec = atciqz(srepr.without_differentials(), astrom)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
srepr = <SphericalRepresentation (lon, lat, distance) in (rad, rad, AU)
(3.50527501, -0.15299223, 0.99699593)>
astrom = (20.78911513, [0.91479605, 0.35636384, 0.15457437], [0.92367677, 0.35156317, 0.15239608], 0.99699593, [-3.97414141e-05...5, 1.41592477e+190, 2.19341577e-152, 1.10639145e+200, 2.34783525e+251, 1.46921998e+195, 3.4430453e-85, 2.65186294e-258)
def atciqz(srepr, astrom):
"""
A slightly modified version of the ERFA function ``eraAtciqz``.
``eraAtciqz`` performs the transformations between two coordinate systems,
with the details of the transformation being encoded into the ``astrom`` array.
There are two issues with the version of atciqz in ERFA. Both are associated
with the handling of light deflection.
The companion function ``eraAticq`` is meant to be its inverse. However, this
is not true for directions close to the Solar centre, since the light deflection
calculations are numerically unstable and therefore not reversible.
This version sidesteps that problem by artificially reducing the light deflection
for directions which are within 90 arcseconds of the Sun's position. This is the
same approach used by the ERFA functions above, except that they use a threshold of
9 arcseconds.
In addition, ERFA's atciqz assumes a distant source, so there is no difference between
the object-Sun vector and the observer-Sun vector. This can lead to errors of up to a
few arcseconds in the worst case (e.g a Venus transit).
Parameters
----------
srepr : `~astropy.coordinates.SphericalRepresentation`
Astrometric ICRS position of object from observer
astrom : eraASTROM array
ERFA astrometry context, as produced by, e.g. ``eraApci13`` or ``eraApcs13``
Returns
--------
ri : float or `~numpy.ndarray`
Right Ascension in radians
di : float or `~numpy.ndarray`
Declination in radians
"""
# ignore parallax effects if no distance, or far away
srepr_distance = srepr.distance
ignore_distance = srepr_distance.unit == u.one
# BCRS coordinate direction (unit vector).
pco = erfa.s2c(srepr.lon.radian, srepr.lat.radian)
# Find BCRS direction of Sun to object
if ignore_distance:
# No distance to object, assume a long way away
q = pco
else:
# Find BCRS direction of Sun to object.
# astrom['eh'] and astrom['em'] contain Sun to observer unit vector,
# and distance, respectively.
eh = astrom['em'][..., np.newaxis] * astrom['eh']
# unit vector from Sun to object
q = eh + srepr_distance[..., np.newaxis].to_value(u.au) * pco
sundist = np.linalg.norm(q, axis=-1, keepdims=True)
> q /= sundist
E RuntimeWarning: invalid value encountered in true_divide
../../.tox/py38-devdeps/lib/python3.8/site-packages/astropy/coordinates/builtin_frames/utils.py:308: RuntimeWarningEven for Sun center, sundist will usually calculate to some very small, but nonzero, distance. I presume that in this freak circumstance, sundist calculated to exactly zero, although I don't understand why NumPy is raising an "invalid value" warning instead of a "divide by zero" warning, so maybe something different has gone wrong?
The obvious workaround is to wrap the problematic division statement with with np.errstate(divide='ignore', invalid='ignore'):.
ping @StuartLittlefair
Activity
StuartLittlefair commentedon Oct 16, 2020
Damn! I don't read my own code sometimes. I even handle the case where sundist is small below. I agree, wrapping with an error state context manager is a good solution, since the
np.wherestatement on the next line will replace any NaNs with the infinite-distance approx.