Repository navigation
Implement direct ICRS -> Observed transform #10887
Description
Description
After years of digging around in the coordinates code I had somehow never realised that we do not implement a direct ICRS -> ALTAZ transform, despite this being perhaps the most commonly used transformation "in the wild". Instead, we apply ICRS -> CIRS -> ALTAZ. This is the way it is done internally inside erfa of course, but there is overhead in doing it this way from the creation of SkyCoord objects and the like.
We should implement the direct transformation within erfa for a speedup of around 6x.
Additional context
A crude attempt at doing this is the following snippet:
import erfa
import numpy as np
from astropy import units as u
from astropy.coordinates.representation import UnitSphericalRepresentation
from astropy.coordinates.transformations import FunctionTransformWithFiniteDifference
from astropy.coordinates.builtin_frames import ICRS, AltAz
from astropy.coordinates.builtin_frames.utils import (get_jd12,
get_polar_motion,
get_cip,
prepare_earth_position_vel)
from astropy.coordinates.baseframe import frame_transform_graph
@frame_transform_graph.transform(FunctionTransformWithFiniteDifference, ICRS, AltAz)
def icrs_to_altaz(icrs_coo, altaz_frame):
jd1_tt, jd2_tt = get_jd12(altaz_frame.obstime, 'tt')
xp, yp = get_polar_motion(altaz_frame.obstime)
x, y, s = get_cip(jd1_tt, jd2_tt)
era = erfa.era00(*get_jd12(altaz_frame.obstime, 'ut1'))
earth_pv, earth_heliocentric = prepare_earth_position_vel(altaz_frame.obstime)
# refraction constants
refa, refb = erfa.refco(
altaz_frame.pressure.to_value(u.hPa),
altaz_frame.temperature.to_value(u.deg_C),
# Relative humidity can be a quantity or a number.
u.Quantity(altaz_frame.relative_humidity).value,
altaz_frame.obswl.to_value(u.micron)
)
astrom = erfa.apco(
jd1_tt, jd2_tt, earth_pv, earth_heliocentric, x, y, s, era,
altaz_frame.location.lon.to_value(u.radian),
altaz_frame.location.lat.to_value(u.radian),
altaz_frame.location.height.to_value(u.m),
xp, yp, s, refa, refb
)
# ignore pm, rv for now, though atciq will handle it for us
ira, idec = icrs_coo.ra.radian, icrs_coo.dec.radian
plx = icrs_coo.distance.to_value(u.arcsecond, equivalencies=u.parallax())
cra, cdec = erfa.atciq(ira, idec, 0, 0, plx, 0, astrom)
az, zen, ha, odec, ora = erfa.atioq(cra, cdec, astrom)
alt = np.pi/2-zen
usrepr = UnitSphericalRepresentation(az*u.radian, alt*u.radian)
return altaz_frame.realize_frame(usrepr)Using this in the following transformation code yields a speedup of 6x on my laptop (from 994ms to 164ms):
import astropy.coordinates as coord
from astropy.time import Time
times = Time('2017-04-06 00:00') + np.arange(0, 1440, 1)*u.min
loc = coord.EarthLocation.of_site('lapalma')
sun = coord.get_body('sun', times, loc)
sun = sun.icrs
aa_frame = coord.AltAz(location=loc, obstime=times)
%time sun.transform_to(aa_frame)Note that a proper implementation should play nicely with the interpolator scheme of @maxnoe, which shouldn't be too hard. And it should play nicely with whatever we do about propagating motion as well.
Activity
maxnoe commentedon Oct 17, 2020
I think this plays into what we discussed in another PR, to move these
get_cipandprepare_..functions intoErfaAstromand always use these. The Interpolator would then provide these higher level functions as compared to the lower-level erfa access.This would make it play nicer with the alternative astrom providers.
StuartLittlefair commentedon Oct 17, 2020
I think it's nice to keep the tools that
erfa_astromuses out oferfa_astrom. But I agree that low-level erfa access should be pruned out of the transformation code and moved intoerfa_astrom.maxnoe commentedon Oct 17, 2020
The point of this interpolation is to make lookup of things faster that are expensive to compute, not to wrap the low-level erfa routines.
It just happened (not by accident, sure) that these were the expensive operations.
But I'd rather have the default provider / interpolating provider have methods like
get_ciporget_earth_pvthan just wrapping the low-level erfa api, since that makes it harder for people implementing the transformations.So the lowest level should probably be the ephemeris provider and upon that we have what we do with the provided values (calculate them for each time, interpolate, do something else).
mkbrewer commentedon Oct 17, 2020
This is exactly what I do in my code. Regarding interpolation, what I do there is to interpolate the output CIRS Ra, Dec coordinates before calling
aper13()to update the earth rotation angle and thenatioq(). That works fine.mkbrewer commentedon Oct 17, 2020
Just a reminder. If you are doing this, you should also have a transformation
icrs_to_hadec(). Either that or justicrs_to_observed()with bothAltAzandHADecattributes since that code yields both.mkbrewer commentedon Oct 18, 2020
@StuartLittlefair I just noticed that there is a bug in your
icrs_to_altaz()function. The third to last argument forapco()should be the TIO locator:Not the CIO locator
s.mkbrewer commentedon Oct 18, 2020
Oh, and the first two arguments for
apco()should betdbrather thantt. I assume that you know that and are intentionally ignoring the ~2 ms difference though. This is entirely reasonable as the time is only used inapcs()to calculate the proper motion intervalastrom['pmt']and stars don't move very far in 2 ms.StuartLittlefair commentedon Oct 19, 2020
Thanks for that @mkbrewer! We are using
ttfor a bit of efficiency since the ERFA docs say it makes no difference for accuracy.I note here for myself when implementing this, that using
erfa.atciqis not a good idea either, because of the issues with light deflection fixed in #10666