Add LVLH frame and geoid intersection method

[jamesgrimmett]

Hi @brandon-rhodes, I made these modifications for a project that I was working on with skyfield, and thought that they might be more widely useful for others in the community.

Changes include;

• Adding an LVLH (Local Vertical Local Horizontal) inertial frame, skyfield.framelib.LVLH, including methods to generate attitude vectors from viewing geometry.
• Adding a method to find the point of intersection between a vector and an ellipsoid, skyfield.toposlib.Geoid.intersection_of

It seems there are a few different conventions for defining an LVLH frame, but I've gone with z-axis anti-aligned with the position vector, and y-axis anti-aligned with the orbital momentum vector.

For the intersection method I've borrowed heavily from here. In short, a "distortion map" is defined to transform the ellipsoid into a unit sphere, and applied to both ellipsoid and vector. Then, the intersection problem is much easier to solve, and the inverse distortion can be applied to the solution.

I'd be interested to hear with you think! If you are in principle keen to add this functionality but have other ideas for the implementation, I'd be interested to discuss that also.

[brandon-rhodes]

I just wanted to leave a quick note to say that I'm planning to look over this by maybe mid-week this week, and respond with some feedback! I'll be interested to see how the code looks.

[jamesgrimmett]

Thanks! No hurry, by the looks of the tests, I still have some tidying up to do.
Also, I thought it might be better if I move 4ca0a99 into a separate PR, so I'll aim to get to that today or tomorrow.

[jamesgrimmett]

Hi @brandon-rhodes, just thought I would ping you again for your thoughts/feedback on this PR. No hurry at all; I'm keen to hear what you think whenever you have the time.

[brandon-rhodes]

Ah, is your tidying-up now complete? Then I'll schedule some time this week to take a look—I'll try to take a look tomorrow morning.

[jamesgrimmett]

Ah yes, sorry about that, I should have dropped you a note to let you know I have tidied up and ran the tests locally to ensure they are passing.

[Ly0n]

This is really a helpful pull request. Thank you so much. I implemented my own LVLH reference frame in the past into skyfield and will compare this results with your implementation. Some feedback here from my side:

1. I would not call it "Looking Vector" but "Line-of-Sight" (LOS).
2. I'm working in Atmospheric Science with limb sounder. Such instruments need a yaw, pitch and roll rotation. The current implantation just has pitch and roll.

In about 2 months I will have a little bit more time working on this. Then I will also be able to open source my implementation to show more practical applications.

[jamesgrimmett]

This is really a helpful pull request. Thank you so much. I implemented my own LVLH reference frame in the past into skyfield and will compare this results with your implementation. Some feedback here from my side:

1. I would not call it "Looking Vector" but "Line-of-Sight" (LOS).

2. I'm working in Atmospheric Science with limb sounder. Such instruments need a yaw, pitch and roll rotation. The current implantation just has pitch and roll.

In about 2 months I will have a little bit more time working on this. Then I will also be able to open source my implementation to show more practical applications.

Thanks! I'm really glad to hear that you have found it useful.

You raise two very good points. I agree with both suggestions, thanks for pointing that out. I was fortunate to only need pitch and roll for the instrument I was working with, so I took a shortcut there. But, completely forgot to go back and add yaw.

I would be very keen to see your implementation/results once they are ready.

I also need to apologise to @brandon-rhodes for making several commits since marking this as ready for review - sorry! I hope it hasn't inconvenienced you. I will aim to make the two changes mentioned above in the next few days, add will add a comment to let you know once I have done that.

[brandon-rhodes]

I also need to apologise to @brandon-rhodes for making several commits since marking this as ready for review - sorry! I hope it hasn't inconvenienced you. I will aim to make the two changes mentioned above in the next few days, add will add a comment to let you know once I have done that.

On the contrary, I need to apologize, for having gotten busy back in June and never coming back to give feedback! I'll try to take a look today, so that you can make adjustments while you're in the code to look at the suggested name improvements.

[brandon-rhodes]

(I did look at the code later that day, by the way, but have not yet grabbed a pad and paper to work through exactly why osculating_elements_of() is being called, and whether there might be a simpler and faster way in Skyfield to get the same vectors. Hopefully later this week!)

[jamesgrimmett]

(I did look at the code later that day, by the way, but have not yet grabbed a pad and paper to work through exactly why osculating_elements_of() is being called, and whether there might be a simpler and faster way in Skyfield to get the same vectors. Hopefully later this week!)

Simpler and faster sounds appealing! I guess somehow we could just use the position and velocity vectors to calculate the frame rotation? Anyway, I won't get ahead of myself, keen to see what you have whenever you have the time.

[brandon-rhodes]

In case I don’t get a chance to go over the math today or tomorrow, I can at least share a few early thoughts about the API rather than making you wait!

• This frame is not inertial — it’s not a non-rotating frame centered at the Solar System Barycenter — so it will be best if it doesn’t inherit from InertialFrame.

• Skyfield has never yet (I don’t think?) had users pass around raw vectors that aren’t part of a bigger object that labels them with units, a center, and a frame of reference. So the final argument to intersection_of() feels a bit raw? I wonder if there’s a way an object could mediate here. Skyfield does try to be minimal, but also to provide context to vectors, like units and orientation.

• It almost feels like we need an object here that represents a vector direction from a specified xyz center. That single object could be returned by a satellite when it’s asked about a given pitch and roll; and then that single object could be consumed by intersection_of().

• Note that in Python, if a method doesn’t use self, then it’s in fact a happy little function and can be removed from its class.

• Let’s put underscores in front of everything, like local_looking_vector(), that we’re not yet directly advertising to users, so we have the leeway to further refine those names in the future without having made promises to user code.

I’ll keep thinking and later this week I’ll hope to make more comments. Thanks for having the PR include tests — that brings it much closer to completion and will make it much easier to merge!

[jamesgrimmett]

It almost feels like we need an object here that represents a vector direction from a specified xyz center. That single object could be returned by a satellite when it’s asked about a given pitch and roll; and then that single object could be consumed by intersection_of()

Perhaps the looking_vector could be refactored as an attitude or line_of_sight function of EarthSatellite, returning an object for intersection_of()? Or perhaps it's better to keep this functionality separate from EarthSatellite? I'll put some thought into this and try put together a more specific description of how this could work. EDIT: scratch that - maybe it's not a bad idea overall, but I think I'd best just focus on representing looking_vector as some kind of vector object for now.

Also, while looking into adding a yaw rotation alongside pitch and roll, I've realised I think I need to be a bit more careful with the way the looking_vector is constructed (i.e., the order that the rotations are applied).

Thanks for the other tips and suggestions - very useful. I will address them and chip away at getting this into a mergeable state.

[jamesgrimmett]

I've had a bit more of a think about your initial feedback, and made some progress (I think..) in the right direction.

• Refactored the functionality from framelib.LVLH.local_looking_vector() to sgp4lib.EarthSatellite._attitude()
  • For now this returns a namedtuple containing the properties of the attitude/line of sight vector. I'm wondering whether a new Attitude(VectorFunction) class would be appropriate here?
• Added a yaw angle alongside pitch and roll, as suggested by @Ly0n
  • I may need to add some more tests to make 100% sure I have implemented this correctly

Interested to hear whether these changes are a step in the right direction, whenever you have some time @brandon-rhodes.

I'll track the progress/remaining to do's in a list here;

• Calculate LVLH frame rotation without calling osculating elements, if possible.
• Add yaw angle to attitude, alongside pitch and roll.
• LVLH frame should not inherit from InertialFrame.
• Define an object to pass to intersection_of.
• Update docstrings and docs after recent refactoring.
• Add underscores in front of everything that we’re not yet directly advertising to users.

[jamesgrimmett]

I've made a set of changes to tick off those remaining items in the to-do list in this comment above. A few points to note/follow up questions;

• It was much simpler than I initially thought to get the rotation matrix of the LVLH frame! I was definitely overthinking that one, no need for the call to osculating_elements_of(), thanks for pointing that out.
• I've used a positionlib.Geometric class (which I now also see is deprecated since I pulled the latest upstream changes) to represent the attitude object returned by EarthSatellite._attitude(). I'm still unsure about whether this is an appropriate choice. Perhaps creating a new Attitude class would be a better choice, or some other existing class that I've overlooked?
• I've added underscores to attitude() and intersection_of() - also just say the word if you'd like me to remove the section that I've added to documentation/earth-satellites.rst if you don't want this functionality to be promised to users yet.

I also rebased this branch to bring in the latest updates to master (i.e., the fixes for the Python 2 tests).

[brandon-rhodes]

@jamesgrimmett — Thank you for the updates! I just this morning have given a conference talk about Skyfield, which will encourage me to now return to the code and look over open issues and pull requests, so you have perfect timing. (It will also help that my spare hours won't be busy with writing the talk!) Hopefully I'll either be able to respond from here on the road, or else next week when I'm home. Feel free to ping me with another comment on the issue if you haven't heard back by the end of next week.

[jamesgrimmett]

Hi @brandon-rhodes - how coincidental! I hope your talk and travel went well.
This morning I was reassuring myself that I had coded the attitude rotation correctly, by comparing results with some existing packages for transforming Euler angles to rotation matrices (e.g., this Matlab package).
Initially I was quite concerned to be getting different results, before I discovered that other packages typically apply intrinsic rotations (reference axes rotate with the body), whereas I have implemented extrinsic rotations (relative to a fixed set of axes). With this is mind, I was then glad to see these results are identical to existing packages. I also learned that applying the same set of rotations intrinsically is as simple as reversing the order of multiplication of the rotation matrices, so that would be simple to implement, if desired.

I have updated the _attitude function docstring to make the extrinsic rotation more explicit, and also removed my usage of the deprecated positionlib.Geometric class. I believe that the testing pipeline should pass since you added a Python 2 fix to jplephem v2.21 (thanks!).

[MichaelBonnet]

If you wish to do some additional applications-oriented testing, I can provide you with some test data generated by a well-known industry tool that includes data like timestamped satellite-Sun LVLH position vectors and associated satellite position/velocity vectors.

[MichaelBonnet]

Instead of or in addition to this method of returning the rotation matrix itself, could it return a scipy Rotation object? it would then be easily convertible back to a rotation matrix or to other useful forms like Euler angles or quaternions.

[jamesgrimmett]

Thanks for reaching out, the test data that you mentioned would be excellent! I had been meaning to spend some time thinking about additional testing strategies, so any data that you could provide would be very much appreciated.
Are you able to attach the data here? Otherwise, could you email it through to - james dot grimmett at protonmail dot com?

I think you make a good point about scipy Rotation objects being a nice, versatile way to carry the rotation information. I would lean toward keeping the return value as a matrix, however. Mainly to remain consistent with several other rotation_at functions scattered throughout the skyfield.framelib module, and also to avoid adding a scipy dependency. That said, I would be curious what @brandon-rhodes opinion is - it may be an option that has been considered before?

[brandon-rhodes]

@jamesgrimmett — I agree entirely. SciPy is even larger than NumPy, and it would be a shame to nearly double the install size of Skyfield just for a single feature by adding it as a dependency. The Rotation object also strikes me as a bit expensive, in terms of its setup cost (as I'm reading through its from_matrix() builder method), for a basic operation in Skyfield.

If there's a good use case for turning Skyfield rotations into Rotation objects, let's just plan on showing in the documentation how it would work, so that folks who need SciPy can import it into their own code and build the Rotation.

[karatemir]

@brandon-rhodes are you planning to merge this? what is the final status?