Add topocentric ITRS frame. Create direct transforms between it and o…

[mkbrewer]

…bserved frames.

Description

This pull request is to address issue #13355. It creates a topocentric ITRS frame and transforms to and from this frame and AltAz or HADec with the ability to add or remove refraction corrections as required.

Closes #13355

Checklist for package maintainer(s)

This checklist is meant to remind the package maintainer(s) who will review this pull request of some common things to look for. This list is not exhaustive.

• Do the proposed changes actually accomplish desired goals?
• Do the proposed changes follow the Astropy coding guidelines?
• Are tests added/updated as required? If so, do they follow the Astropy testing guidelines?
• Are docs added/updated as required? If so, do they follow the Astropy documentation guidelines?
• Is rebase and/or squash necessary? If so, please provide the author with appropriate instructions. Also see "When to rebase and squash commits".
• Did the CI pass? If no, are the failures related? If you need to run daily and weekly cron jobs as part of the PR, please apply the Extra CI label.
• Is a change log needed? If yes, did the change log check pass? If no, add the no-changelog-entry-needed label. If this is a manual backport, use the skip-changelog-checks label unless special changelog handling is necessary.
• Is this a big PR that makes a "What's new?" entry worthwhile and if so, is (1) a "what's new" entry included in this PR and (2) the "whatsnew-needed" label applied?
• Is a milestone set? Milestone must be set but astropy-bot check might be missing; do not let the green checkmark fool you.
• At the time of adding the milestone, if the milestone set requires a backport to release branch(es), apply the appropriate backport-X.Y.x label(s) before merge.

[pep8speaks]

Hello @mkbrewer 👋! It looks like you've made some changes in your pull request, so I've checked the code again for style.

There are no PEP8 style issues with this pull request - thanks! 🎉

Comment last updated at 2022-07-27 19:50:52 UTC

[mkbrewer]

I am concerned that the transform graph in the Read the Docs build under Astronomical Coordinate Systems Built-in Frame Classes doesn't show my new ITRS<->AltAz and ITRS<->HADec transforms. Am I missing something there?

[mkbrewer]

I have been doing some more testing. My direct ITRS<->Observed transforms speed up transforms to AltAz from ITRS, TETE, and TEME by a factor of 4. TETE seems to automatically find the shortest path through the ITRS rather than through GCRS as before. For TEME, I used the example given in Working with Earth Satellites Using Astropy Coordinates in the documentation. That example shows a transformation from TEME to AltAz. The result is incorrect as it naturally assumes that the TEME coordinates include aberration as seen from the geocenter, which near Earth satellites do not. I found the following offsets after doing the transformation correctly as described in my updated documentation for the ITRS frame:

Az:   59.28803392 deg delta: 46.0031arcsec
El:   10.94799670 deg delta: 15.472arcsec
Dist: 1431.96640195 km delta: 4.723 m

I will update the example to show the correct way of doing this. BTW, it would be nice if the units label for arcsec included a space after the number.

[mkbrewer]

Could someone please find the time to review this? @mhvk? @StuartLittlefair? I am hoping that it won't just languish for months or years.

[StuartLittlefair]

Hi @mkbrewer - yes of course, sorry for leaving it to gather dust. I'll try and get to this early next week, but please remind me if I don't!

[StuartLittlefair]

Hi @mkbrewer - got around to this finally.

I like the code, but wonder what we can do to make this more user friendly. As you can see in my inline comments, we are straying into the area where the astropy user really has to understand what the coordinate transforms do in order to use them correctly. I think most of this can be addressed in a separate PR for documentation, but perhaps you could look at some of your new docstrings in this light as well?

[StuartLittlefair]

Hi @mkbrewer

Thanks for the careful replies. I'm definitely in favour of this PR, and I think some minor tweaks to the docs will help the end user. I've suggested these in the inline comments.

There are some minor code style errors that are preventing that CI test from passing. You need 2 blank lines around L292 and L954 in test_intermediate_transformations.py.

Finally, I think this is a change that would benefit from highlighting in the "What's New" section of the 5.2 release. Something along the lines of "A topocentric ITRS frame has been added that makes dealing with near-Earth objects easier" followed by a code example like:

    t = Time('J2010')
    obj = EarthLocation(-1*u.deg, 52*u.deg, height=10.*u.km)
    home = EarthLocation(-1*u.deg, 52*u.deg, height=0.*u.km)

    # Direction of object from GEOCENTER
    itrs_geo = obj.get_itrs(t).cartesian

    # now get the Geocentric ITRS position of observatory
    obsrepr = home.get_itrs(t).cartesian

    # topocentric ITRS position of a straight overhead object
    itrs_repr = itrs_geo - obsrepr

    # create a ITRS object that appears straight overhead for a TOPOCENTRIC OBSERVER
    topocentric_itrs_frame = ITRS(obstime=t, location=home)
    itrs_topo = topocentric_itrs_frame.realize_frame(itrs_repr)

    # convert to AltAz
    aa = itrs_topo.transform_to(AltAz(obstime=t, location=home))

[mkbrewer]

One further question that I have: #13398 (comment)

[StuartLittlefair]

One further question that I have: #13398 (comment)

Very good question! Yes - you should be concerned. I can't believe I missed this, but the transforms don't get registered with the graph unless the code in itrs_observed_transforms.py is executed, which happens at import time. You need to import this module in astropy/coordinates/builtin_frames/__init__.py for your new transforms to be registered!

[mkbrewer]

Thanks!

[mkbrewer]

Still no luck with the transform graph.

I don't know what it is about subtracting the two Cartesian positions in satellite.rst. This used to work, now it complains about differentials.

[StuartLittlefair]

Well, I don't know about the graphic in the docs, but your transform is there in the frame transform graph. Previously this used to go through CIRS:

In [41]: from astropy.coordinates.builtin_frames import frame_transform_graph
    ...: from astropy.coordinates import AltAz, ITRS
    ...: 
    ...: trans = frame_transform_graph.get_transform(ITRS, AltAz)
    ...: for t in trans.transforms:
    ...:     print(f'{t.fromsys.name} -> {t.tosys.name}')
    ...: 
itrs -> altaz

For the line in satellites.rst, to be honest I don't know how this ever managed to work, but you need to replace L108 with

topo_itrs_repr = itrs_geo.cartesian.without_differentials() - siding_spring.get_itrs(t).cartesian

or, if you want to keep the satellite velocity (which you get using transform_to):

itrs_siding_spring = siding_spring.get_gcrs(t).transform_to(ITRS(obstime=t))   # CLUNKY, should add zero velocity to get_itrs
dp = itrs_geo.cartesian.without_differentials() - itrs_siding_spring.cartesian.without_differentials()
dv = itrs_geo.cartesian.differentials['s'] - itrs_siding_spring.cartesian.differentials['s']
topo_itrs_repr = dp.with_differentials(dv)

N.B - perhaps I'm overdoing this, but is it worth explaining in satellites.rst why you shouldn't just use itrs_geo.transform_to(AltAz(obstime=t, location=siding_spring))?

[mkbrewer]

Oh yes. I somehow failed to notice the velocity in the TEME frame and didn't include that in my testing, which is why it worked. Velocity in the ITRS is independent of location, so I don't see the point of dv here. Why not just:

topo_itrs_repr = dp.with_differentials(itrs_geo.cartesian.differentials['s'])

[mkbrewer]

Regarding the graphic. I found this in index.rst:

.. automodapi:: astropy.coordinates.builtin_frames
    :skip: make_transform_graph_docs
    :no-inheritance-diagram:

Does that mean building a new graph is being skipped?

[mkbrewer]

I'm going to go with an example of carrying the velocity to the topocentric frame without dv.

[mkbrewer]

Looks good now. The one failure is unrelated.

[StuartLittlefair]

Lovely! This looks almost perfect now. I have one more small request, which is that you add a test of the round-trip ITRS->HADec->ITRS. This will test L142 in itrs_observed_transforms.py which is currently untested

[mkbrewer]

Unfortunately, the documentation build failed, but whatever caused that appears to be unrelated to this PR. And this after my transforms finally appeared in the transform graph after the last build.

[pllim]

Hmm... WARNING: The following HTTPError has occurred fetching https://matplotlib.org/stable//: 522 () -- I'll restart RTD build and see.

[pllim]

Is the coordinates failure in https://github.com/astropy/astropy/runs/7472108451?check_suite_focus=true related?

The matplotlib stuff looks unrelated.

[mkbrewer]

No, it is not.

[pllim]

Yeah, ok. I see it in main too.

As far as CI is concerned, I think this is good, unless you also want "exotic" archs to run, then we can add "Extra CI" label; lemme know. Thanks!

[StuartLittlefair]

Looks great @mkbrewer. Thanks for the effort; it cleans up dealing with near earth objects substantially. I'll hold off merging for a few days in case @mhvk wants to add anything.

[mhvk]

This looks great! All my comments are essentially nitpicks, so I'm approving, but perhaps this gives the opportunity to also squash the commits (at least the bug/typo fixes).

The new examples make me wonder whether it makes sense to allow subtraction of ITRS frames... it seems a bit ugly to have to go through the cartesian coordinates. But that is another issue!

[mhvk]

Hmm, reminds me that this erfa function needs a quantity-aware wrapper as well! But not yet there, so good as is.

[mkbrewer]

You guys get my taskmaster of the year award! 😜

[mkbrewer]

Failures are unrelated again.

[mhvk]

@mkbrewer - OK, looks all OK to go in. Just squashed the commits a bit, to remove the PEP ones, etc. Will merge once the tests ran again. Thanks so much for a really nice contribution!

[mkbrewer]

You are welcome! I hope that this is helpful for those who deal with the positions of satellites.