Skip to content

Commit 60edc62

Browse files
embrayembray
committed
Merge pull request #4228 from embray/fits/issue-3854
String-related fixes for FITS tables
1 parent 3e82e35 commit 60edc62

10 files changed

Lines changed: 710 additions & 457 deletions

‎CHANGES.rst‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -101,6 +101,11 @@ Bug Fixes
101101

102102
- ``astropy.coordinates``
103103

104+
- Fixed errors in transformations for objects within a few AU of the
105+
Earth. Included substansive changes to transformation machinery
106+
that may change distances at levels ~machine precision for other
107+
objects. [#4254]
108+
104109
- ``astropy.cosmology``
105110

106111
- ``astropy.io.ascii``

‎astropy/coordinates/baseframe.py‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -28,7 +28,6 @@
2828
SphericalRepresentation,
2929
UnitSphericalRepresentation,
3030
REPRESENTATION_CLASSES)
31-
from .earth import EarthLocation
3231

3332

3433
__all__ = ['BaseCoordinateFrame', 'frame_transform_graph', 'GenericFrame',
@@ -1213,3 +1212,7 @@ def __setattr__(self, name, value):
12131212
raise AttributeError("can't set frame attribute '{0}'".format(name))
12141213
else:
12151214
super(GenericFrame, self).__setattr__(name, value)
1215+
1216+
# doing this import at the bottom prevents a circular import issue that is
1217+
# otherwise present due to EarthLocation needing to import ITRS
1218+
from .earth import EarthLocation

‎astropy/coordinates/builtin_frames/cirs_observed_transforms.py‎

Lines changed: 41 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,8 @@
1212
from ... import units as u
1313
from ..baseframe import frame_transform_graph
1414
from ..transformations import FunctionTransform
15-
from ..representation import SphericalRepresentation, UnitSphericalRepresentation
15+
from ..representation import (SphericalRepresentation, CartesianRepresentation,
16+
UnitSphericalRepresentation)
1617
from ... import _erfa as erfa
1718

1819
from .cirs import CIRS
@@ -22,21 +23,25 @@
2223

2324
@frame_transform_graph.transform(FunctionTransform, CIRS, AltAz)
2425
def cirs_to_altaz(cirs_coo, altaz_frame):
25-
if np.all(cirs_coo.obstime != altaz_frame.obstime):
26+
if np.any(cirs_coo.obstime != altaz_frame.obstime):
2627
# the only frame attribute for the current CIRS is the obstime, but this
2728
# would need to be updated if a future change allowed specifying an
2829
# Earth location algorithm or something
2930
cirs_coo = cirs_coo.transform_to(CIRS(obstime=altaz_frame.obstime))
30-
srepr = cirs_coo.represent_as(UnitSphericalRepresentation)
31-
cirs_ra = srepr.lon.to(u.radian).value
32-
cirs_dec = srepr.lat.to(u.radian).value
31+
32+
# we use the same obstime everywhere now that we know they're the same
33+
obstime = cirs_coo.obstime
34+
35+
usrepr = cirs_coo.represent_as(UnitSphericalRepresentation)
36+
cirs_ra = usrepr.lon.to(u.radian).value
37+
cirs_dec = usrepr.lat.to(u.radian).value
3338

3439
lon, lat, height = altaz_frame.location.to_geodetic('WGS84')
35-
xp, yp = get_polar_motion(cirs_coo.obstime)
40+
xp, yp = get_polar_motion(obstime)
3641

37-
#first set up the astrometry context for ICRS<->CIRS
38-
astrom = erfa.apio13(cirs_coo.obstime.jd1, cirs_coo.obstime.jd2,
39-
get_dut1utc(cirs_coo.obstime),
42+
#first set up the astrometry context for CIRS<->AltAz
43+
astrom = erfa.apio13(obstime.jd1, obstime.jd2,
44+
get_dut1utc(obstime),
4045
lon.to(u.radian).value, lat.to(u.radian).value,
4146
height.to(u.m).value,
4247
xp, yp, # polar motion
@@ -54,9 +59,10 @@ def cirs_to_altaz(cirs_coo, altaz_frame):
5459
lon=u.Quantity(az, u.radian, copy=False),
5560
copy=False)
5661
else:
57-
#now we get the distance as just the cartesian
58-
#distance from the earth location to the coordinate location
59-
distance = altaz_frame.location.itrs.separation_3d(cirs_coo)
62+
# now we get the distance as the cartesian distance from the earth
63+
# location to the coordinate location
64+
locitrs = altaz_frame.location.get_itrs(obstime)
65+
distance = locitrs.separation_3d(cirs_coo)
6066
rep = SphericalRepresentation(lat=u.Quantity(PIOVER2 - zen, u.radian, copy=False),
6167
lon=u.Quantity(az, u.radian, copy=False),
6268
distance=distance,
@@ -66,9 +72,9 @@ def cirs_to_altaz(cirs_coo, altaz_frame):
6672

6773
@frame_transform_graph.transform(FunctionTransform, AltAz, CIRS)
6874
def altaz_to_cirs(altaz_coo, cirs_frame):
69-
srepr = altaz_coo.represent_as(UnitSphericalRepresentation)
70-
az = srepr.lon.to(u.radian).value
71-
zen = PIOVER2 - srepr.lat.to(u.radian).value
75+
usrepr = altaz_coo.represent_as(UnitSphericalRepresentation)
76+
az = usrepr.lon.to(u.radian).value
77+
zen = PIOVER2 - usrepr.lat.to(u.radian).value
7278

7379
lon, lat, height = altaz_coo.location.to_geodetic('WGS84')
7480
xp, yp = get_polar_motion(altaz_coo.obstime)
@@ -86,17 +92,32 @@ def altaz_to_cirs(altaz_coo, cirs_frame):
8692
altaz_coo.obswl.value)
8793

8894
# the 'A' indicates zen/az inputs
89-
cirs_ra, cirs_dec = erfa.atoiq('A', az, zen, astrom)
95+
cirs_ra, cirs_dec = erfa.atoiq('A', az, zen, astrom)*u.radian
9096

9197
if isinstance(altaz_coo.data, UnitSphericalRepresentation):
9298
distance = None
9399
else:
94-
#now we get the distance as just the cartesian
95-
#distance from the earth location to the coordinate location
96-
distance = altaz_coo.separation_3d(altaz_coo.location.itrs)
100+
locitrs = altaz_coo.location.get_itrs(altaz_coo.obstime)
101+
102+
# To compute the distance in a way that is reversable with cirs_to_altaz
103+
# we use basic trigonometry. The altaz_coo's distance ("d") is one leg
104+
# of a triangle, and the earth center -> EarthLocation distance ("r")
105+
# is a neighboring leg. We can also easily get the angle between the
106+
# earth center->target (calculated above by apio13), and
107+
# earth center->EarthLocation vectors. This is a Side-Side-Angle
108+
# situation, and the formula below is the trig formula to solve for
109+
# the remaining side
110+
111+
ucirs = cirs_frame.realize_frame(UnitSphericalRepresentation(lon=cirs_ra, lat=cirs_dec))
112+
delta = ucirs.separation(locitrs)
113+
r = locitrs.spherical.distance
114+
d = altaz_coo.distance
115+
116+
sindoverd = np.sin(delta) / d
117+
distance = np.sin(delta + np.arcsin(r*sindoverd))/sindoverd
97118

98119
#the final transform may be a no-op if the obstimes are the same
99-
return CIRS(ra=cirs_ra*u.radian, dec=cirs_dec*u.radian, distance=distance,
120+
return CIRS(ra=cirs_ra, dec=cirs_dec, distance=distance,
100121
obstime=altaz_coo.obstime).transform_to(cirs_frame)
101122

102123

‎astropy/coordinates/builtin_frames/icrs_cirs_transforms.py‎

Lines changed: 101 additions & 81 deletions
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,8 @@
1212
from ... import units as u
1313
from ..baseframe import frame_transform_graph
1414
from ..transformations import FunctionTransform
15-
from ..representation import UnitSphericalRepresentation, SphericalRepresentation
15+
from ..representation import (SphericalRepresentation, CartesianRepresentation,
16+
UnitSphericalRepresentation)
1617
from ... import _erfa as erfa
1718

1819
from .icrs import ICRS
@@ -24,61 +25,70 @@
2425

2526
@frame_transform_graph.transform(FunctionTransform, ICRS, CIRS)
2627
def icrs_to_cirs(icrs_coo, cirs_frame):
27-
#parallax in arcsec
28-
if isinstance(icrs_coo.data, UnitSphericalRepresentation): # no distance
29-
px = 0
30-
else:
31-
px = 1 / icrs_coo.distance.to(u.parsec).value
32-
srepr = icrs_coo.represent_as(UnitSphericalRepresentation)
33-
i_ra = srepr.lon.to(u.radian).value
34-
i_dec = srepr.lat.to(u.radian).value
35-
3628
#first set up the astrometry context for ICRS<->CIRS
3729
astrom, eo = erfa.apci13(cirs_frame.obstime.jd1, cirs_frame.obstime.jd2)
3830

39-
# TODO: possibly switch to something that is like atciq, but skips the first
40-
# step, b/c that involves some wasteful computations b/c pm=0 here
41-
cirs_ra, cirs_dec = erfa.atciq(i_ra, i_dec, 0, 0, px, 0, astrom)
31+
if icrs_coo.data.get_name() == 'unitspherical' or icrs_coo.data.to_cartesian().x.unit == u.one:
32+
# if no distance, just do the infinite-distance/no parallax calculation
33+
usrepr = icrs_coo.represent_as(UnitSphericalRepresentation)
34+
i_ra = usrepr.lon.to(u.radian).value
35+
i_dec = usrepr.lat.to(u.radian).value
36+
cirs_ra, cirs_dec = erfa.atciqz(i_ra, i_dec, astrom)
4237

43-
dat = icrs_coo.data
44-
if dat.get_name() == 'unitspherical' or dat.to_cartesian().x.unit == u.one:
45-
rep = UnitSphericalRepresentation(lat=u.Quantity(cirs_dec, u.radian, copy=False),
46-
lon=u.Quantity(cirs_ra, u.radian, copy=False),
47-
copy=False)
38+
newrep = UnitSphericalRepresentation(lat=u.Quantity(cirs_dec, u.radian, copy=False),
39+
lon=u.Quantity(cirs_ra, u.radian, copy=False),
40+
copy=False)
4841
else:
49-
#compute the distance as just the cartesian sum from moving to the Earth
50-
#we have to do this because the ERFA functions throw away distance info
51-
distance = np.sum((astrom['eb']*u.au + icrs_coo.cartesian.xyz.T)**2, -1)**0.5
52-
rep = SphericalRepresentation(lat=u.Quantity(cirs_dec, u.radian, copy=False),
53-
lon=u.Quantity(cirs_ra, u.radian, copy=False),
54-
distance=distance, copy=False)
42+
# When there is a distance, we first offset for parallax to get the
43+
# astrometric coordinate direction and *then* run the ERFA transform for
44+
# no parallax/PM. This ensures reversiblity and is more sensible for
45+
# inside solar system objects
46+
newxyz = icrs_coo.cartesian.xyz.T - astrom['eb']*u.au
47+
newcart = CartesianRepresentation(newxyz.T)
48+
49+
srepr = newcart.represent_as(SphericalRepresentation)
50+
i_ra = srepr.lon.to(u.radian).value
51+
i_dec = srepr.lat.to(u.radian).value
52+
cirs_ra, cirs_dec = erfa.atciqz(i_ra, i_dec, astrom)
5553

56-
return cirs_frame.realize_frame(rep)
54+
newrep = SphericalRepresentation(lat=u.Quantity(cirs_dec, u.radian, copy=False),
55+
lon=u.Quantity(cirs_ra, u.radian, copy=False),
56+
distance=srepr.distance, copy=False)
57+
58+
return cirs_frame.realize_frame(newrep)
5759

5860
@frame_transform_graph.transform(FunctionTransform, CIRS, ICRS)
5961
def cirs_to_icrs(cirs_coo, icrs_frame):
6062
srepr = cirs_coo.represent_as(UnitSphericalRepresentation)
6163
cirs_ra = srepr.lon.to(u.radian).value
6264
cirs_dec = srepr.lat.to(u.radian).value
6365

64-
#first set up the astrometry context for ICRS<->CIRS
66+
# set up the astrometry context for ICRS<->cirs and then convert to
67+
# astrometric coordinate direction
6568
astrom, eo = erfa.apci13(cirs_coo.obstime.jd1, cirs_coo.obstime.jd2)
69+
i_ra, i_dec = erfa.aticq(cirs_ra, cirs_dec, astrom)
70+
71+
if cirs_coo.data.get_name() == 'unitspherical' or cirs_coo.data.to_cartesian().x.unit == u.one:
72+
# if no distance, just use the coordinate direction to yield the
73+
# infinite-distance/no parallax answer
74+
newrep = UnitSphericalRepresentation(lat=u.Quantity(i_dec, u.radian, copy=False),
75+
lon=u.Quantity(i_ra, u.radian, copy=False),
76+
copy=False)
77+
else:
78+
# When there is a distance, apply the parallax/offset to the SSB as the
79+
# last step - ensures round-tripping with the icrs_to_cirs transform
6680

67-
icrs_ra, icrs_dec = erfa.aticq(cirs_ra, cirs_dec, astrom)
81+
# the distance in intermedrep is *not* a real distance as it does not
82+
# include the offset back to the SSB
83+
intermedrep = SphericalRepresentation(lat=u.Quantity(i_dec, u.radian, copy=False),
84+
lon=u.Quantity(i_ra, u.radian, copy=False),
85+
distance=cirs_coo.distance,
86+
copy=False)
6887

69-
dat = cirs_coo.data
70-
if dat.get_name() == 'unitspherical' or dat.to_cartesian().x.unit == u.one:
71-
rep = UnitSphericalRepresentation(lat=u.Quantity(icrs_dec, u.radian, copy=False),
72-
lon=u.Quantity(icrs_ra, u.radian, copy=False),
73-
copy=False)
74-
else:
75-
#compute the distance as just the cartesian sum from moving to the SSB
76-
#we have to do this because the ERFA functions throw away distance info
77-
distance = np.sum((-astrom['eb']*u.au + cirs_coo.cartesian.xyz.T)**2, -1)**0.5
78-
rep = SphericalRepresentation(lat=u.Quantity(icrs_dec, u.radian, copy=False),
79-
lon=u.Quantity(icrs_ra, u.radian, copy=False),
80-
distance=distance, copy=False)
81-
return icrs_frame.realize_frame(rep)
88+
newxyz = intermedrep.to_cartesian().xyz.T + astrom['eb']*u.au
89+
newrep = CartesianRepresentation(newxyz.T).represent_as(SphericalRepresentation)
90+
91+
return icrs_frame.realize_frame(newrep)
8292

8393
@frame_transform_graph.transform(FunctionTransform, CIRS, CIRS)
8494
def cirs_to_cirs(from_coo, to_frame):
@@ -98,65 +108,75 @@ def cirs_to_cirs(from_coo, to_frame):
98108

99109
@frame_transform_graph.transform(FunctionTransform, ICRS, GCRS)
100110
def icrs_to_gcrs(icrs_coo, gcrs_frame):
101-
#parallax in arcsec
102-
if isinstance(icrs_coo.data, UnitSphericalRepresentation): # no distance
103-
px = 0
104-
else:
105-
px = 1 / icrs_coo.distance.to(u.parsec).value
106-
srepr = icrs_coo.represent_as(UnitSphericalRepresentation)
107-
i_ra = srepr.lon.to(u.radian).value
108-
i_dec = srepr.lat.to(u.radian).value
109-
110111
#first set up the astrometry context for ICRS<->GCRS
111112
pv = np.array([gcrs_frame.obsgeoloc.value,
112113
gcrs_frame.obsgeovel.value])
113114
astrom = erfa.apcs13(gcrs_frame.obstime.jd1, gcrs_frame.obstime.jd2, pv)
114115

115-
# TODO: possibly switch to something that is like atciq, but skips the first
116-
# step, b/c that involves some wasteful computations b/c pm=0 here
117-
gcrs_ra, gcrs_dec = erfa.atciq(i_ra, i_dec, 0, 0, px, 0, astrom)
116+
if icrs_coo.data.get_name() == 'unitspherical' or icrs_coo.data.to_cartesian().x.unit == u.one:
117+
# if no distance, just do the infinite-distance/no parallax calculation
118+
usrepr = icrs_coo.represent_as(UnitSphericalRepresentation)
119+
i_ra = usrepr.lon.to(u.radian).value
120+
i_dec = usrepr.lat.to(u.radian).value
121+
gcrs_ra, gcrs_dec = erfa.atciqz(i_ra, i_dec, astrom)
118122

119-
dat = icrs_coo.data
120-
if dat.get_name() == 'unitspherical' or dat.to_cartesian().x.unit == u.one:
121-
rep = UnitSphericalRepresentation(lat=u.Quantity(gcrs_dec, u.radian, copy=False),
122-
lon=u.Quantity(gcrs_ra, u.radian, copy=False),
123-
copy=False)
123+
newrep = UnitSphericalRepresentation(lat=u.Quantity(gcrs_dec, u.radian, copy=False),
124+
lon=u.Quantity(gcrs_ra, u.radian, copy=False),
125+
copy=False)
124126
else:
125-
#compute the distance as just the cartesian sum from moving to the Earth
126-
#we have to do this because the ERFA functions throw away distance info
127-
distance = np.sum((astrom['eb']*u.au + icrs_coo.cartesian.xyz.T)**2, -1)**0.5
128-
rep = SphericalRepresentation(lat=u.Quantity(gcrs_dec, u.radian, copy=False),
129-
lon=u.Quantity(gcrs_ra, u.radian, copy=False),
130-
distance=distance, copy=False)
131-
return gcrs_frame.realize_frame(rep)
127+
# When there is a distance, we first offset for parallax to get the
128+
# BCRS coordinate direction and *then* run the ERFA transform for no
129+
# parallax/PM. This ensures reversiblity and is more sensible for
130+
# inside solar system objects
131+
newxyz = icrs_coo.cartesian.xyz.T - astrom['eb']*u.au
132+
newcart = CartesianRepresentation(newxyz.T)
133+
134+
srepr = newcart.represent_as(SphericalRepresentation)
135+
i_ra = srepr.lon.to(u.radian).value
136+
i_dec = srepr.lat.to(u.radian).value
137+
gcrs_ra, gcrs_dec = erfa.atciqz(i_ra, i_dec, astrom)
138+
139+
newrep = SphericalRepresentation(lat=u.Quantity(gcrs_dec, u.radian, copy=False),
140+
lon=u.Quantity(gcrs_ra, u.radian, copy=False),
141+
distance=srepr.distance, copy=False)
142+
143+
return gcrs_frame.realize_frame(newrep)
132144

133145

134146
@frame_transform_graph.transform(FunctionTransform, GCRS, ICRS)
135147
def gcrs_to_icrs(gcrs_coo, icrs_frame):
136148
srepr = gcrs_coo.represent_as(UnitSphericalRepresentation)
137-
cirs_ra = srepr.lon.to(u.radian).value
138-
cirs_dec = srepr.lat.to(u.radian).value
149+
gcrs_ra = srepr.lon.to(u.radian).value
150+
gcrs_dec = srepr.lat.to(u.radian).value
139151

140-
#first set up the astrometry context for ICRS<->GCRS
152+
# set up the astrometry context for ICRS<->GCRS and then convert to BCRS
153+
# coordinate direction
141154
pv = np.array([gcrs_coo.obsgeoloc.value,
142155
gcrs_coo.obsgeovel.value])
143156
astrom = erfa.apcs13(gcrs_coo.obstime.jd1, gcrs_coo.obstime.jd2, pv)
157+
i_ra, i_dec = erfa.aticq(gcrs_ra, gcrs_dec, astrom)
158+
159+
if gcrs_coo.data.get_name() == 'unitspherical' or gcrs_coo.data.to_cartesian().x.unit == u.one:
160+
# if no distance, just use the coordinate direction to yield the
161+
# infinite-distance/no parallax answer
162+
newrep = UnitSphericalRepresentation(lat=u.Quantity(i_dec, u.radian, copy=False),
163+
lon=u.Quantity(i_ra, u.radian, copy=False),
164+
copy=False)
165+
else:
166+
# When there is a distance, apply the parallax/offset to the SSB as the
167+
# last step - ensures round-tripping with the icrs_to_gcrs transform
144168

145-
icrs_ra, icrs_dec = erfa.aticq(cirs_ra, cirs_dec, astrom)
169+
# the distance in intermedrep is *not* a real distance as it does not
170+
# include the offset back to the SSB
171+
intermedrep = SphericalRepresentation(lat=u.Quantity(i_dec, u.radian, copy=False),
172+
lon=u.Quantity(i_ra, u.radian, copy=False),
173+
distance=gcrs_coo.distance,
174+
copy=False)
146175

147-
dat = gcrs_coo.data
148-
if dat.get_name() == 'unitspherical' or dat.to_cartesian().x.unit == u.one:
149-
rep = UnitSphericalRepresentation(lat=u.Quantity(icrs_dec, u.radian, copy=False),
150-
lon=u.Quantity(icrs_ra, u.radian, copy=False),
151-
copy=False)
152-
else:
153-
#compute the distance as just the cartesian sum from moving to the SSB
154-
#we have to do this because the ERFA functions throw away distance info
155-
distance = np.sum((astrom['eb']*u.au + gcrs_coo.cartesian.xyz.T)**2, axis=-1)**0.5
156-
rep = SphericalRepresentation(lat=u.Quantity(icrs_dec, u.radian, copy=False),
157-
lon=u.Quantity(icrs_ra, u.radian, copy=False),
158-
distance=distance, copy=False)
159-
return icrs_frame.realize_frame(rep)
176+
newxyz = intermedrep.to_cartesian().xyz.T + astrom['eb']*u.au
177+
newrep = CartesianRepresentation(newxyz.T).represent_as(SphericalRepresentation)
178+
179+
return icrs_frame.realize_frame(newrep)
160180

161181

162182
@frame_transform_graph.transform(FunctionTransform, GCRS, GCRS)

‎astropy/coordinates/builtin_frames/itrs.py‎

Lines changed: 0 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -16,9 +16,6 @@ class ITRS(BaseCoordinateFrame):
1616
defined by a series of reference locations near the surface of the Earth.
1717
For more background on the ITRS, see the references provided in the
1818
:ref:`astropy-coordinates-seealso` section of the documentation.
19-
20-
Note that this frame does *not* have an ``obstime``, because it is
21-
rotating with the Earth.
2219
"""
2320

2421
default_representation = CartesianRepresentation

0 commit comments

Comments
 (0)