Returning GCRS Coordinates for Moon. Solves poliastro#382

src/poliastro/twobody/orbit.py

@@ -5,7 +5,7 @@
 from astropy import units as u
 from astropy import time
 
-from astropy.coordinates import CartesianRepresentation, get_body_barycentric_posvel, get_body_barycentric
+from astropy.coordinates import CartesianRepresentation, get_body_barycentric_posvel, get_body_barycentric, ICRS, GCRS, CartesianDifferential, solar_system_ephemeris
 
 from poliastro.constants import J2000
 from poliastro.twobody.angles import nu_to_M, E_to_nu
@@ -15,6 +15,8 @@
 from poliastro.twobody import classical
 from poliastro.twobody import equinoctial
 
+from poliastro.bodies import Moon
+
 from ._base import BaseState
 
 
@@ -162,7 +164,24 @@ def from_body_ephem(cls, body, epoch=None):
             warn("Input time was converted to scale='tdb' with value "
                  "{}. Use Time(..., scale='tdb') instead."
                  .format(epoch.tdb.value), TimeScaleWarning)
+        if body == Moon:
+            solar_system_ephemeris.set("jpl")
         r, v = get_body_barycentric_posvel(body.name, epoch)
+        if body == Moon:
+            moon_icrs = ICRS(
+                x=r.x, y=r.y, z=r.z,
+                v_x=v.x, v_y=v.y, v_z=v.z,
+                representation=CartesianRepresentation,
+                differential_type=CartesianDifferential
+            )
+            moon_gcrs = moon_icrs.transform_to(GCRS(obstime=epoch))
+            moon_gcrs.representation = CartesianRepresentation
+            return cls.from_vectors(
+                body.parent,
+                [moon_gcrs.x, moon_gcrs.y, moon_gcrs.z] * u.km,
+                [moon_gcrs.v_x, moon_gcrs.v_y, moon_gcrs.v_z] * (u.km / u.s),
+                epoch=epoch
+            )
         return cls.from_vectors(body.parent, r.xyz.to(u.km), v.xyz.to(u.km / u.day), epoch)
 
     @classmethod