Add `BaseBodycentricRepresentation` #14851

cmarmo · 2023-05-19T21:07:01Z

Description

This pull request adds the possibility to create geodetic representations described by [planetocentric latitudes]
~~(https://en.wikipedia.org/wiki/Latitude#Geodetic_and_geocentric_latitudes), west positive longitudes, longitudes spanning from 0 to 360 degrees.~~
The default behavior is backward compatible.
Suggested in #11170 (comment).

Towards #11170.

This work is funded by the Europlanet 2024 Research Infrastructure (RI) Grant.

github-actions · 2023-05-19T21:07:34Z

cmarmo · 2023-05-19T22:34:36Z

It seems like I'm still confused about units conversions... :(

cmarmo · 2023-05-25T20:14:10Z

For the sake of completeness, this pull request makes the BaseGeodeticRepresentation closer to the Planetary Data System class Geodetic_Model.

mhvk

@cmarmo - Good to have these options, but I do think at least some of them should be subclasses rather than options on the base class. For _ellipsoid', I ended up agreeing it was not necessary, since having it on the baseclass gives a little overhead only at class creation time, and in the end it is really the same thing as far as the actual calculations are concerned. But here the calculations for conversion to/from cartesian are effectively separated, making the code more complicated and giving overhead for every calculation. Generally, I think each different interpretation of the coordinates should have its own (base) class -- just like PhysicsSpherical is a different class from Spherical. I would say the rule should be that the moment from_cartesian and to_cartesian depend on any class attribute, it would be better to have a different class.
So, I think we definitely should have a separate base class for the planetocentric latitudes (BasePlanetoCentric?). I think we also will be better off with a separate class for West longitudes, though given the rarity of its use, perhaps this is best implemented via a mixin class that overrides from_cartesian and to_cartesian changing the sign of the longitude appropriately.
Only for _wrap_angle, I think it is fine to leave that on the base class -- the calculations do not depend on it, one can see it as similar to changing the class of the attributes.

astropy/coordinates/representation/geodetic.py

mhvk · 2023-05-26T12:16:59Z

astropy/coordinates/representation/geodetic.py

@@ -45,9 +56,16 @@ class BaseGeodeticRepresentation(BaseRepresentation):
    to quantities holding correct values (with units of length and dimensionless,
    respectively), or alternatively an ``_ellipsoid`` attribute to the relevant ERFA
    index (as passed in to `erfa.eform`).
+    Longitudes are east positive and span from -180 to 180 degrees by default.


I think we should have separate classes; see main comment.

astropy/coordinates/representation/geodetic.py

astropy/coordinates/tests/test_geodetic_representations.py

astropy/coordinates/tests/test_representation.py

astropy/coordinates/representation/__init__.py

mhvk

Thanks for making the change to a different base class. That does feel much better!

But, seeing this again, I wonder about the approach for wrap_angle. Really, the problem is that we have a Longitude class that has a default of 360 degrees, while we would like 180 degrees. Would it make sense to just define

class Longitude180(Longitude):
    _default_wrap_angle = Angle(180 * u.deg)

and then use that? Then anyone who wants to have a different default, can just choose to use a different Longitude class. And we don't have to deal with the new _wrap_angle attribute.

p.s. Need to think a bit about the mixin still... And I haven't looked in detail over the tests.

astropy/coordinates/representation/geodetic.py

mhvk · 2023-06-14T16:16:06Z

astropy/coordinates/representation/geodetic.py

@@ -62,12 +98,21 @@ def __init_subclass__(cls, **kwargs):
            raise AttributeError(
                f"{cls.__name__} requires '_ellipsoid' or '_equatorial_radius' and '_flattening'."
            )
+        if not hasattr(


Small detail, but could do

if not u.Quantity(cls._wrap_angle).unit.is_equivalent(u.deg):

(but see my larger comment about _wrap_angle)

astropy/coordinates/representation/geodetic.py

astropy/coordinates/tests/test_representation.py

cmarmo · 2023-06-15T08:52:13Z

Thanks for your review @mhvk !
I have addressed the "obvious" fixes.
I have to better understand units and wrap angles for the rest of the comments.... I will be back ASAP.
Thanks again for your precious comments.

cmarmo · 2023-06-16T22:33:56Z

The tests are failing because the class WGS84BodycentricRepresentation, which I have defined in order to test the ability of using _ellipsoid with a BaseBodycentricRepresentation, persists and it is used by earth.py.
My understanding was that the teardown_function would take care of "preserving the original REPRESENTATION_CLASSES", so this failure should not happen.
Anyway it also seems like the selection of a representation only based on the ellipsoid in EarthLocation may lead to mistakes now that BaseBodycentricRepresentation also uses an ellipsoid.

On the long term, I see the importance of changing EarthLocation to a representation, but how should I proceed on the short term? Should I remove the possibility of defining a bodycentric representation using an ellipsoid?

Any suggestion? Thanks a lot!

mhvk · 2023-06-17T20:10:31Z

Simplest right now would seem to not allow ellipsoid for Bodycentric. As is, ellipsoid is only defined for Earth, and it is only used in EarthLocation, so I think this is fine.

cmarmo · 2023-06-20T01:04:32Z

Hi @mhvk . I believe I have removed all the unnecessary conversions to degrees, following your comments.
I have updated the title of the pull request and the changelog to reflect what is done in the code: ie the introduction of a new BaseBodycentricRepresentation.
However, I have kept the attribute _wrap_angle and defaulted it at 360 for BaseBodycentricRepresentation as this is the default recommended by Cartographic Coordinates & Rotational Elements Working Group (see https://rdcu.be/b32WL page 26)....
I'm not completely convinced in having it removed ...

But if we stick with using _wrap_angle, we need to think a bit about how we can potentially make it faster - right now, the wrapping is done twice (I think).

I guess this is also because I don't understand this previous comment

mhvk · 2023-06-22T17:35:38Z

@cmarmo - thanks for the changes. I'm still not sold on the _wrap_angle attribute - essentially you are making how data is displayed part of the class, while I feel what the class should cares about is only how data is used. And for the use of a longitude, it really doesn't matter what the wrap angle is (just like its units don't matter). I still think the way forward is to define a new Longitude180 class and use that where relevant.

However, I also think this discussion is actually distracting from the main point of this PR, which is to introduce a BaseBodycentricRepresentation class. How about splitting this up in pieces? I.e., here remove all stuff with wrap_angle and just get the basic class working.

Indeed, while I think the WestLongitude mixin class is the right approach, I'd even tend to put that in a separate PR as well, where it could include an update to the docs that make clearer what the use case is.

cmarmo · 2023-06-22T20:33:21Z

However, I also think this discussion is actually distracting from the main point of this PR, which is to introduce a BaseBodycentricRepresentation class. How about splitting this up in pieces? I.e., here remove all stuff with wrap_angle and just get the basic class working.

Indeed, while I think the WestLongitude mixin class is the right approach, I'd even tend to put that in a separate PR as well, where it could include an update to the docs that make clearer what the use case is.

I guess I'm loosing focus indeed... I need the Bodycentric representation for the frame definition, while west/east and 180/360 are not needed.... and I'm running out of time... ⏲️
Let's rescope the PR.

mhvk · 2023-06-22T21:09:46Z

OK, thanks, ping me when I should look again. And apologies for not realizing this earlier!

cmarmo · 2023-06-22T21:23:26Z

And apologies for not realizing this earlier!

No problem! I'm treasuring all the comments for the follow-up! :)

cmarmo · 2023-06-22T22:35:08Z

@mhvk thanks for your patience!
Is the PR acceptable in its current status?
Should I add something in the what's new section?

mhvk

This looks great! There are only nitpicks and suggestions for the docstrings left. But also: do indeed add to the whatsnew entry you already created for other geodetic representations.

mhvk · 2023-06-22T23:09:53Z

astropy/coordinates/representation/geodetic.py

        lon, lat, height = erfa.gc2gde(
            cls._equatorial_radius, cls._flattening, cart.get_xyz(xyz_axis=-1)
        )
-        return cls(lon, lat, height, copy=False)
+        return cls(


Nitpick, but could you get this back to the default format? Be sure to delete the final , after copy=False, otherwise black will break it over several lines again.

mhvk · 2023-06-22T23:10:09Z

astropy/coordinates/representation/geodetic.py

@@ -45,6 +46,7 @@ class BaseGeodeticRepresentation(BaseRepresentation):
    to quantities holding correct values (with units of length and dimensionless,
    respectively), or alternatively an ``_ellipsoid`` attribute to the relevant ERFA
    index (as passed in to `erfa.eform`).
+    Longitudes are east positive and span from 0 to 360 degrees by default.


Nitpick: might as well leave this out for now.

Instead, I think here is the place to mention how lon, lat, and height are defined relative to the ellipsoid.

Alternatively, you could clarify the definition in geodetic_base_doc above (and then not use it for bodycentric)

mhvk · 2023-06-22T23:12:41Z

astropy/coordinates/representation/geodetic.py

+    Subclasses need to set attributes ``_equatorial_radius`` and ``_flattening``
+    to quantities holding correct values (with units of length and dimensionless,
+    respectively).
+    Longitudes are east positive and span from 0 to 360 degrees by default.


Again, I'd remove this for now and instead describe how the parameters are used (and how it is different from geodetic).

Alternatively, do not use geodetic_base_doc, but instead write a complete docstring here that also describes how lat is defined.

mhvk · 2023-06-22T23:14:52Z

astropy/coordinates/representation/geodetic.py

+            lon,
+            lat,
+            height,
+            copy=False,


Again, could you bring this on a single line?

mhvk · 2023-06-22T23:49:00Z

docs/changes/coordinates/14851.feature.rst

@@ -0,0 +1,2 @@
+Add ``BaseBodycentricRepresentation``, a new spheroidal representation for bodycentric
+latitudes and longitudes spanning from 0 to 360 degrees.


Let's remove the "spanning from 0 to 360 degrees" here too.

mhvk · 2023-06-22T23:50:28Z

docs/coordinates/representations.rst

+`~astropy.coordinates.BaseBodycentricRepresentation` is the coordinate representation
+recommended by the Cartographic Coordinates & Rotational Elements Working Group
+(see for example its `2019 report <https://rdcu.be/b32WL>`_): the bodicentric latitude
+and longitude are spherical latitude and longitude originated at the barycenter of the


originated at -> relative to

docs/coordinates/representations.rst

mhvk · 2023-06-22T23:53:20Z

docs/coordinates/representations.rst

@@ -48,6 +50,13 @@ between the vertical to the surface at a specific point of the spheroid and its
 projection onto the equatorial plane.
 The latitude is a value ranging from -90 to 90 degrees, the longitude from 0 to 360
 degrees.
+`~astropy.coordinates.BaseBodycentricRepresentation` is the coordinate representation
+recommended by the Cartographic Coordinates & Rotational Elements Working Group
+(see for example its `2019 report <https://rdcu.be/b32WL>`_): the bodicentric latitude


bodicentric -> bodycentric

docs/coordinates/representations.rst

eerovaher

I'm letting @mhvk do the main review, but I do have a few comments.

eerovaher · 2023-06-23T15:45:14Z

astropy/coordinates/representation/geodetic.py

+        x_spheroid = self._equatorial_radius * coslat * coslon
+        y_spheroid = self._equatorial_radius * coslat * sinlon
+        z_spheroid = self._equatorial_radius * (1 - self._flattening) * sinlat


These variables are unused.

Suggested change

x_spheroid = self._equatorial_radius * coslat * coslon

y_spheroid = self._equatorial_radius * coslat * sinlon

z_spheroid = self._equatorial_radius * (1 - self._flattening) * sinlat

Indeed! They are no longer needed as the calculation is done for r already.

eerovaher · 2023-06-23T15:50:04Z

astropy/coordinates/representation/geodetic.py

+        z_spheroid = self._equatorial_radius * (1 - self._flattening) * sinlat
+        r = (
+            self._equatorial_radius
+            * np.sqrt(coslat**2 + ((1 - self._flattening) * sinlat) ** 2)


Could use np.hypot() here:

Suggested change

* np.sqrt(coslat**2 + ((1 - self._flattening) * sinlat) ** 2)

* np.hypot(coslat, (1 - self._flattening) * sinlat)

eerovaher · 2023-06-23T15:57:05Z

astropy/coordinates/tests/test_geodetic_representations.py

+    [
+        WGS84GeodeticRepresentation,
+        IAUMARS2000BodycentricRepresentation,
+    ],


Doesn't this fit on one line?

Suggested change

[

WGS84GeodeticRepresentation,

IAUMARS2000BodycentricRepresentation,

],

[WGS84GeodeticRepresentation, IAUMARS2000BodycentricRepresentation],

If it does then the other similar cases should be updated analogously.

eerovaher · 2023-06-23T16:00:40Z

docs/coordinates/representations.rst

+degrees, the height is the elevation above the surface of the spheroid (measured on the
+perpendicular to the surface).


Should this be replaced with the following?

Suggested change

degrees, the height is the elevation above the surface of the spheroid (measured on the

perpendicular to the surface).

degrees, the height is the elevation above the surface of the spheroid (measured

perpendicular to the surface).

mhvk

@cmarmo - only some nitpicks left beyond @eerovaher's comments. With those, I'll be happy to merge!

mhvk · 2023-06-23T17:18:10Z

astropy/coordinates/representation/geodetic.py

+
+    Subclasses need to set attributes ``_equatorial_radius`` and ``_flattening``
+    to quantities holding correct values (with units of length and dimensionless,
+    respectively). the bodicentric latitude and longitude are spherical latitude


"the bodicentric" -> "The bodycentric"

mhvk · 2023-06-23T17:20:02Z

astropy/coordinates/representation/geodetic.py

+        x_spheroid = self._equatorial_radius * coslat * coslon
+        y_spheroid = self._equatorial_radius * coslat * sinlon
+        z_spheroid = self._equatorial_radius * (1 - self._flattening) * sinlat


Indeed! They are no longer needed as the calculation is done for r already.

mhvk · 2023-06-23T17:21:27Z

astropy/coordinates/representation/geodetic.py

+        d = np.hypot(p, cart.z)
+        lat = np.arctan2(cart.z, p)
+        p_spheroid = cls._equatorial_radius * np.cos(lat)
+        z_spheroid = (cls._equatorial_radius * (1 - cls._flattening)) * np.sin(lat)


Total nitpick, but outer parentheses are not necessary.

mhvk · 2023-06-23T17:24:55Z

docs/coordinates/representations.rst

 create specific representations on spheroidal bodies.
-This is used internally for the standard Earth ellipsoids used in
+`~astropy.coordinates.BaseGeodeticRepresentation` is used internally for the standard
+Earth ellipsoids used in
 `~astropy.coordinates.EarthLocation`


Maybe combine with last line while you are making the last changes.

Sorry I'm not sure to understand this comment....

Meant to combine the two lines to give```
Earth ellipsoids used in ~astropy.coordinates.EarthLocation

mhvk · 2023-06-23T17:28:04Z

docs/coordinates/representations.rst

+recommended by the Cartographic Coordinates & Rotational Elements Working Group
+(see for example its `2019 report <https://rdcu.be/b32WL>`_): the bodycentric latitude
+and longitude are spherical latitude and longitude relative to the barycenter of the
+body, the height is the distance from the spheroid surface on the line of sight.


"the height is the distance from the spheroid surface on the line of sight" -> "height is the distance from the surface (measured radially)."

mhvk · 2023-06-23T17:28:31Z

docs/coordinates/representations.rst

@@ -706,11 +717,13 @@ In pseudo-code, this means that a class will look like::

 .. _astropy-coordinates-create-geodetic:

-Creating Your Own Geodetic Representation
-----------------------------------------
+Creating Your Own Geodetic and Bodycentric Representation


Maybe make plural: "Representations"

mhvk · 2023-06-23T17:46:49Z

@cmarmo - separately from this PR, my sense would be to try to follow this up with adding some of the I/O and WCS system things you have been mentioning rather than worry about west longitude and wrap angles!

mhvk · 2023-06-23T19:54:25Z

docs/coordinates/representations.rst

 create specific representations on spheroidal bodies.
-This is used internally for the standard Earth ellipsoids used in


I didn't mean to have this part about the internal use to be removed!

mhvk · 2023-06-23T19:55:32Z

astropy/coordinates/representation/geodetic.py

@@ -110,7 +110,7 @@ class BaseBodycentricRepresentation(BaseRepresentation):

    Subclasses need to set attributes ``_equatorial_radius`` and ``_flattening``
    to quantities holding correct values (with units of length and dimensionless,
-    respectively). the bodicentric latitude and longitude are spherical latitude
+    respectively). the bodycentric latitude and longitude are spherical latitude


Still need to capitalize "The"

mhvk · 2023-06-23T19:56:36Z

docs/coordinates/representations.rst

 The latitude is a value ranging from -90 to 90 degrees, the longitude from 0 to 360
 degrees.

+`~astropy.coordinates.BaseGeodeticRepresentation` is used internally for the standard


Oops, sorry, you moved it here! This is better!

mhvk

This all looks great now. I think I'll just merge, the one uncapitalized "t" can be done at some other time!

mhvk · 2023-06-23T19:59:42Z

Thanks, @cmarmo!!

cmarmo · 2023-06-23T20:13:18Z

Thanks @mhvk for your patience!

github-actions bot added the coordinates label May 19, 2023

cmarmo mentioned this pull request May 19, 2023

Planetary WCS implementation in Astropy #11170

Closed

pllim added this to the v6.0 milestone May 19, 2023

pllim requested a review from mhvk May 19, 2023 22:02

cmarmo force-pushed the geodetic-repr-options branch from d3f53fd to a6e30fb Compare May 19, 2023 22:26

cmarmo force-pushed the geodetic-repr-options branch 3 times, most recently from f1a8a58 to 25fb355 Compare May 25, 2023 20:11

mhvk requested changes May 26, 2023

View reviewed changes

cmarmo force-pushed the geodetic-repr-options branch from 25fb355 to 50c29f4 Compare June 14, 2023 10:50

cmarmo commented Jun 14, 2023

View reviewed changes

astropy/coordinates/representation/__init__.py Show resolved Hide resolved

mhvk reviewed Jun 14, 2023

View reviewed changes

cmarmo changed the title ~~Add planetary useful options to BaseGeodeticRepresentation~~ Add BaseBodycentricRepresentation Jun 16, 2023

cmarmo force-pushed the geodetic-repr-options branch from d17d2dc to 265e9bb Compare June 19, 2023 21:09

Add bodycentric representation. Add tests and documentation.

918a6b3

cmarmo force-pushed the geodetic-repr-options branch from 317bc2b to 918a6b3 Compare June 22, 2023 21:53

mhvk reviewed Jun 22, 2023

View reviewed changes

Address comments.

8eda98e

eerovaher reviewed Jun 23, 2023

View reviewed changes

mhvk reviewed Jun 23, 2023

View reviewed changes

Address comments.

6a5ad91

mhvk reviewed Jun 23, 2023

View reviewed changes

mhvk approved these changes Jun 23, 2023

View reviewed changes

mhvk added the Enhancement label Jun 23, 2023

mhvk merged commit 9e3d836 into astropy:main Jun 23, 2023
29 checks passed

cmarmo deleted the geodetic-repr-options branch June 23, 2023 20:13

		@@ -0,0 +1,2 @@
		Add ``BaseBodycentricRepresentation``, a new spheroidal representation for bodycentric
		latitudes and longitudes spanning from 0 to 360 degrees.

	x_spheroid = self._equatorial_radius * coslat * coslon
	y_spheroid = self._equatorial_radius * coslat * sinlon
	z_spheroid = self._equatorial_radius * (1 - self._flattening) * sinlat

	* np.sqrt(coslat*2 + ((1 - self._flattening) sinlat) ** 2)
	* np.hypot(coslat, (1 - self._flattening) * sinlat)

		degrees, the height is the elevation above the surface of the spheroid (measured on the
		perpendicular to the surface).

		create specific representations on spheroidal bodies.
		This is used internally for the standard Earth ellipsoids used in

Add BaseBodycentricRepresentation #14851

Add BaseBodycentricRepresentation #14851

Conversation

cmarmo commented May 19, 2023 • edited

Description

github-actions bot commented May 19, 2023 • edited by mhvk

cmarmo commented May 19, 2023

cmarmo commented May 25, 2023

mhvk left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mhvk left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

cmarmo commented Jun 15, 2023

cmarmo commented Jun 16, 2023

mhvk commented Jun 17, 2023

cmarmo commented Jun 20, 2023

mhvk commented Jun 22, 2023

cmarmo commented Jun 22, 2023

mhvk commented Jun 22, 2023

cmarmo commented Jun 22, 2023

cmarmo commented Jun 22, 2023

mhvk left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

eerovaher left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mhvk left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mhvk commented Jun 23, 2023

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mhvk left a comment

Choose a reason for hiding this comment

mhvk commented Jun 23, 2023

cmarmo commented Jun 23, 2023

Add `BaseBodycentricRepresentation` #14851

Add `BaseBodycentricRepresentation` #14851

cmarmo commented May 19, 2023 •

edited

github-actions bot commented May 19, 2023 •

edited by mhvk