Skip to content

Propagating objects in time as part of the transformation framework #10404

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
ayshih wants to merge 4 commits into from
Closed

Propagating objects in time as part of the transformation framework #10404

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
11c9d08
First draft of propagators for ICRS<->BCRS transformations
ayshih May 26, 2020
10e8031
Fix for scalar data and array obstime
ayshih May 26, 2020
b7ba882
Miscellaneous fixes
ayshih May 27, 2020
7708e03
Make BasePropagator a proper abstract base class
ayshih May 28, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions astropy/coordinates/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,3 +23,4 @@
from .solar_system import *
from .spectral_quantity import *
from .spectral_coordinate import *
from .propagators import *
3 changes: 2 additions & 1 deletion astropy/coordinates/builtin_frames/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@
from .cirs import CIRS
from .itrs import ITRS
from .hcrs import HCRS
from .bcrs import BCRS
from .teme import TEME
from .ecliptic import * # there are a lot of these so we don't list them all explicitly
from .skyoffset import SkyOffsetFrame
Expand All @@ -59,7 +60,7 @@
# get included
__all__ = ['ICRS', 'FK5', 'FK4', 'FK4NoETerms', 'Galactic', 'Galactocentric',
'galactocentric_frame_defaults',
'Supergalactic', 'AltAz', 'GCRS', 'CIRS', 'ITRS', 'HCRS', 'TEME',
'Supergalactic', 'AltAz', 'GCRS', 'CIRS', 'ITRS', 'HCRS', 'BCRS', 'TEME',
'PrecessedGeocentric', 'GeocentricMeanEcliptic',
'BarycentricMeanEcliptic', 'HeliocentricMeanEcliptic',
'GeocentricTrueEcliptic', 'BarycentricTrueEcliptic',
Expand Down
29 changes: 29 additions & 0 deletions astropy/coordinates/builtin_frames/bcrs.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,29 @@
# -*- coding: utf-8 -*-
# Licensed under a 3-clause BSD style license - see LICENSE.rst

from astropy.utils.decorators import format_doc
from astropy.coordinates.attributes import TimeAttribute
from .utils import DEFAULT_OBSTIME
from astropy.coordinates.baseframe import base_doc
from .baseradec import BaseRADecFrame, doc_components

__all__ = ['BCRS']


doc_footer = """
Other parameters
----------------
obstime : `~astropy.time.Time`
The time at which the observation is taken.
"""


@format_doc(base_doc, components=doc_components, footer=doc_footer)
class BCRS(BaseRADecFrame):
"""
Like ICRS, but with an ``obstime`` attribute
"""

obstime = TimeAttribute(default=DEFAULT_OBSTIME)

# Transformations are defined in icrs_cirs_transforms.py
2 changes: 1 addition & 1 deletion astropy/coordinates/builtin_frames/hcrs.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,4 +42,4 @@ class HCRS(BaseRADecFrame):

obstime = TimeAttribute(default=DEFAULT_OBSTIME)

# Transformations are defined in icrs_circ_transforms.py
# Transformations are defined in icrs_cirs_transforms.py
90 changes: 56 additions & 34 deletions astropy/coordinates/builtin_frames/icrs_cirs_transforms.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,9 @@

from astropy import units as u
from astropy.coordinates.baseframe import frame_transform_graph
from astropy.coordinates.transformations import FunctionTransformWithFiniteDifference, AffineTransform
from astropy.coordinates.propagators import transform_propagator
from astropy.coordinates.transformations import (FunctionTransformWithFiniteDifference,
FunctionTransform, AffineTransform)
from astropy.coordinates.representation import (SphericalRepresentation,
CartesianRepresentation,
UnitSphericalRepresentation,
Expand All @@ -20,22 +22,42 @@
from .gcrs import GCRS
from .cirs import CIRS
from .hcrs import HCRS
from .bcrs import BCRS
from .utils import get_jd12, aticq, atciqz, get_cip, prepare_earth_position_vel


# First the ICRS/CIRS related transforms
@frame_transform_graph.transform(FunctionTransformWithFiniteDifference, ICRS, CIRS)
def icrs_to_cirs(icrs_coo, cirs_frame):
# first set up the astrometry context for ICRS<->CIRS
@frame_transform_graph.transform(FunctionTransform, ICRS, BCRS)
def icrs_to_bcrs(icrs_coo, bcrs_frame):
if transform_propagator.get():
return transform_propagator.get().propagate(icrs_coo, bcrs_frame.obstime)
else:
data = icrs_coo.data
if data.isscalar and not bcrs_frame.obstime.isscalar:

@mhvk mhvk May 27, 2020

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yikes! This probably is the best argument yet for allowing attributes that have a bigger shape than the data, and then just broadcasting the data - but almost certainly best not in this PR! See #8819 for an earlier issue about this. I think we should fix this in the initializer...

data = data._apply('repeat', bcrs_frame.obstime.shape)
return bcrs_frame.realize_frame(data)


@frame_transform_graph.transform(FunctionTransform, BCRS, ICRS)
def bcrs_to_icrs(bcrs_coo, icrs_frame):
if transform_propagator.get():
return transform_propagator.get().depropagate(bcrs_coo)
else:
return icrs_frame.realize_frame(bcrs_coo.data)


# First the BCRS/CIRS related transforms
@frame_transform_graph.transform(FunctionTransformWithFiniteDifference, BCRS, CIRS)
def bcrs_to_cirs(bcrs_coo, cirs_frame):
# first set up the astrometry context for BCRS<->CIRS
jd1, jd2 = get_jd12(cirs_frame.obstime, 'tt')
x, y, s = get_cip(jd1, jd2)
earth_pv, earth_heliocentric = prepare_earth_position_vel(cirs_frame.obstime)
# erfa.apci requests TDB but TT can be used instead of TDB without any significant impact on accuracy
astrom = erfa.apci(jd1, jd2, earth_pv, earth_heliocentric, x, y, s)

if icrs_coo.data.get_name() == 'unitspherical' or icrs_coo.data.to_cartesian().x.unit == u.one:
if bcrs_coo.data.get_name() == 'unitspherical' or bcrs_coo.data.to_cartesian().x.unit == u.one:
# if no distance, just do the infinite-distance/no parallax calculation
usrepr = icrs_coo.represent_as(UnitSphericalRepresentation)
usrepr = bcrs_coo.represent_as(UnitSphericalRepresentation)
i_ra = usrepr.lon.to_value(u.radian)
i_dec = usrepr.lat.to_value(u.radian)
cirs_ra, cirs_dec = atciqz(i_ra, i_dec, astrom)
Expand All @@ -50,7 +72,7 @@ def icrs_to_cirs(icrs_coo, cirs_frame):
# inside solar system objects
astrom_eb = CartesianRepresentation(astrom['eb'], unit=u.au,
xyz_axis=-1, copy=False)
newcart = icrs_coo.cartesian - astrom_eb
newcart = bcrs_coo.cartesian - astrom_eb

srepr = newcart.represent_as(SphericalRepresentation)
i_ra = srepr.lon.to_value(u.radian)
Expand All @@ -64,13 +86,13 @@ def icrs_to_cirs(icrs_coo, cirs_frame):
return cirs_frame.realize_frame(newrep)


@frame_transform_graph.transform(FunctionTransformWithFiniteDifference, CIRS, ICRS)
def cirs_to_icrs(cirs_coo, icrs_frame):
@frame_transform_graph.transform(FunctionTransformWithFiniteDifference, CIRS, BCRS)
def cirs_to_bcrs(cirs_coo, bcrs_frame):
srepr = cirs_coo.represent_as(SphericalRepresentation)
cirs_ra = srepr.lon.to_value(u.radian)
cirs_dec = srepr.lat.to_value(u.radian)

# set up the astrometry context for ICRS<->cirs and then convert to
# set up the astrometry context for BCRS<->cirs and then convert to
# astrometric coordinate direction
jd1, jd2 = get_jd12(cirs_coo.obstime, 'tt')
x, y, s = get_cip(jd1, jd2)
Expand All @@ -87,7 +109,7 @@ def cirs_to_icrs(cirs_coo, icrs_frame):
copy=False)
else:
# When there is a distance, apply the parallax/offset to the SSB as the
# last step - ensures round-tripping with the icrs_to_cirs transform
# last step - ensures round-tripping with the bcrs_to_cirs transform

# the distance in intermedrep is *not* a real distance as it does not
# include the offset back to the SSB
Expand All @@ -100,7 +122,7 @@ def cirs_to_icrs(cirs_coo, icrs_frame):
xyz_axis=-1, copy=False)
newrep = intermedrep + astrom_eb

return icrs_frame.realize_frame(newrep)
return bcrs_frame.realize_frame(newrep)


@frame_transform_graph.transform(FunctionTransformWithFiniteDifference, CIRS, CIRS)
Expand All @@ -117,11 +139,11 @@ def cirs_to_cirs(from_coo, to_frame):
return from_coo.transform_to(ICRS).transform_to(to_frame)


# Now the GCRS-related transforms to/from ICRS
# Now the GCRS-related transforms to/from BCRS

@frame_transform_graph.transform(FunctionTransformWithFiniteDifference, ICRS, GCRS)
def icrs_to_gcrs(icrs_coo, gcrs_frame):
# first set up the astrometry context for ICRS<->GCRS. There are a few steps...
@frame_transform_graph.transform(FunctionTransformWithFiniteDifference, BCRS, GCRS)
def bcrs_to_gcrs(bcrs_coo, gcrs_frame):
# first set up the astrometry context for BCRS<->GCRS. There are a few steps...
# get the position and velocity arrays for the observatory. Need to
# have xyz in last dimension, and pos/vel in one-but-last.
# (Note could use np.stack once our minimum numpy version is >=1.10.)
Expand All @@ -136,9 +158,9 @@ def icrs_to_gcrs(icrs_coo, gcrs_frame):
# get astrometry context object, astrom.
astrom = erfa.apcs(jd1, jd2, obs_pv, earth_pv, earth_heliocentric)

if icrs_coo.data.get_name() == 'unitspherical' or icrs_coo.data.to_cartesian().x.unit == u.one:
if bcrs_coo.data.get_name() == 'unitspherical' or bcrs_coo.data.to_cartesian().x.unit == u.one:
# if no distance, just do the infinite-distance/no parallax calculation
usrepr = icrs_coo.represent_as(UnitSphericalRepresentation)
usrepr = bcrs_coo.represent_as(UnitSphericalRepresentation)
i_ra = usrepr.lon.to_value(u.radian)
i_dec = usrepr.lat.to_value(u.radian)
gcrs_ra, gcrs_dec = atciqz(i_ra, i_dec, astrom)
Expand All @@ -153,7 +175,7 @@ def icrs_to_gcrs(icrs_coo, gcrs_frame):
# inside solar system objects
astrom_eb = CartesianRepresentation(astrom['eb'], unit=u.au,
xyz_axis=-1, copy=False)
newcart = icrs_coo.cartesian - astrom_eb
newcart = bcrs_coo.cartesian - astrom_eb

srepr = newcart.represent_as(SphericalRepresentation)
i_ra = srepr.lon.to_value(u.radian)
Expand All @@ -168,13 +190,13 @@ def icrs_to_gcrs(icrs_coo, gcrs_frame):


@frame_transform_graph.transform(FunctionTransformWithFiniteDifference,
GCRS, ICRS)
def gcrs_to_icrs(gcrs_coo, icrs_frame):
GCRS, BCRS)
def gcrs_to_bcrs(gcrs_coo, bcrs_frame):
srepr = gcrs_coo.represent_as(SphericalRepresentation)
gcrs_ra = srepr.lon.to_value(u.radian)
gcrs_dec = srepr.lat.to_value(u.radian)

# set up the astrometry context for ICRS<->GCRS and then convert to BCRS
# set up the astrometry context for BCRS<->GCRS and then convert to BCRS
# coordinate direction
obs_pv = erfa.pav2pv(
gcrs_coo.obsgeoloc.get_xyz(xyz_axis=-1).to_value(u.m),
Expand All @@ -195,7 +217,7 @@ def gcrs_to_icrs(gcrs_coo, icrs_frame):
copy=False)
else:
# When there is a distance, apply the parallax/offset to the SSB as the
# last step - ensures round-tripping with the icrs_to_gcrs transform
# last step - ensures round-tripping with the bcrs_to_gcrs transform

# the distance in intermedrep is *not* a real distance as it does not
# include the offset back to the SSB
Expand All @@ -208,7 +230,7 @@ def gcrs_to_icrs(gcrs_coo, icrs_frame):
xyz_axis=-1, copy=False)
newrep = intermedrep + astrom_eb

return icrs_frame.realize_frame(newrep)
return bcrs_frame.realize_frame(newrep)


@frame_transform_graph.transform(FunctionTransformWithFiniteDifference, GCRS, GCRS)
Expand All @@ -235,7 +257,7 @@ def gcrs_to_hcrs(gcrs_coo, hcrs_frame):
gcrs_ra = srepr.lon.to_value(u.radian)
gcrs_dec = srepr.lat.to_value(u.radian)

# set up the astrometry context for ICRS<->GCRS and then convert to ICRS
# set up the astrometry context for BCRS<->GCRS and then convert to BCRS
# coordinate direction
obs_pv = erfa.pav2pv(
gcrs_coo.obsgeoloc.get_xyz(xyz_axis=-1).to_value(u.m),
Expand Down Expand Up @@ -280,12 +302,12 @@ def gcrs_to_hcrs(gcrs_coo, hcrs_frame):

_NEED_ORIGIN_HINT = ("The input {0} coordinates do not have length units. This "
"probably means you created coordinates with lat/lon but "
"no distance. Heliocentric<->ICRS transforms cannot "
"no distance. Heliocentric<->BCRS transforms cannot "
"function in this case because there is an origin shift.")


@frame_transform_graph.transform(AffineTransform, HCRS, ICRS)
def hcrs_to_icrs(hcrs_coo, icrs_frame):
@frame_transform_graph.transform(AffineTransform, HCRS, BCRS)
def hcrs_to_bcrs(hcrs_coo, bcrs_frame):
# this is just an origin translation so without a distance it cannot go ahead
if isinstance(hcrs_coo.data, UnitSphericalRepresentation):
raise u.UnitsError(_NEED_ORIGIN_HINT.format(hcrs_coo.__class__.__name__))
Expand All @@ -305,13 +327,13 @@ def hcrs_to_icrs(hcrs_coo, icrs_frame):
return None, bary_sun_pos


@frame_transform_graph.transform(AffineTransform, ICRS, HCRS)
def icrs_to_hcrs(icrs_coo, hcrs_frame):
@frame_transform_graph.transform(AffineTransform, BCRS, HCRS)
def bcrs_to_hcrs(bcrs_coo, hcrs_frame):
# this is just an origin translation so without a distance it cannot go ahead
if isinstance(icrs_coo.data, UnitSphericalRepresentation):
raise u.UnitsError(_NEED_ORIGIN_HINT.format(icrs_coo.__class__.__name__))
if isinstance(bcrs_coo.data, UnitSphericalRepresentation):
raise u.UnitsError(_NEED_ORIGIN_HINT.format(bcrs_coo.__class__.__name__))

if icrs_coo.data.differentials:
if bcrs_coo.data.differentials:
from astropy.coordinates.solar_system import get_body_barycentric_posvel
bary_sun_pos, bary_sun_vel = get_body_barycentric_posvel('sun',
hcrs_frame.obstime)
Expand Down
Loading