You have most of the pieces here already.
UPDATED
The proper coordinates
In this problem, we want to get the positions of astronomical objects in terms of a Cartesian system that is geocentric and rotates with the Earth. In astronomy, the positions of objects are commonly given in terms of right ascension (RA) and declination (Dec), which are similar to longitude and latitude projected onto the sky, but do not rotate with the Earth. Because of this similarity RA and Dec are natural coordinates to start with: it is simple to take the RA and Dec and then boost them into a frame that rotates with the Earth. One can then make a simple transformation from spherical into Cartesian coordinates.
$$
(r, Dec, RA) \rightarrow (r, lat, long) \rightarrow (x, y, x)
$$
Aside:
Because the Earth's axis moves slowly with respect the distant stars (precession and nutation etc.), there is some ambiguity as to whether RA and Dec should move with the Earth or remain fixed. As detailed at length in this question, the various astronomical Data functions when applied to the Sun, the Moon, the planets, and the other planetary moons supply precessing (true geocentric) right ascension and declination. But, for stars besides the Sun, StarData provides a non-precessing (fixed) RA and Dec. All appear to be accurate up to about a degree or RA and less than a degree of Dec. I haven't checked yet what the other astronomical Data functions do, so I will restrict this answer to the Sun, the Moon, and planets. It should be sufficiently easy to generalize.
Conversion into Earth-rotating coordinates
To convert from geocentric, precessing RA into a coordinate in an Earth-rotating reference frame, we calculate the local hour angle as
$$
LHA(long,t) = LST(long,t) - RA(t)
$$
Here, $LST$ is the local sidereal time. Since we are interested in the position of the Sun or Moon with respect to Earth longitude and latitude, we pick our earth position to be $(lat, long) = (0, 0)$. Conveniently, the local sidereal time, on the prime meridian is Greenwich sidereal time (GST, closely related to GMT), so a "longitude-type" angle (only "longitude-type", because we don't need to require this angle to only run between -180 degrees and 180 degrees) for the moon or sun is given as
$$
LHA(0,t) = 360^\circ - \left( GST(t) - RA(t) \right)
$$
(the 360 degrees converts the angle from westward-directed to eastward-directed, so that we can get a right-handed coordinate system, with the x-axis sticking out from the equator-prime meridian point, the y-axis going off through the equator at 90 degrees, and the z-axis along up the axis of rotation)
Mathematica implementation
We can put this into Mathematica now fairly easily:
Sun position
SunPositionXYZ[date_] := With[
{date0 = TimeZoneConvert[date, 0], pos0 = GeoPosition[{0, 0}]},
Module[{pos, GST, r, lat, long},
pos = StarData["Sun",
{EntityProperty["Star", "RightAscension", {"Date" -> date0}],
EntityProperty["Star", "Declination", {"Date" -> date0}]}];
pos = UnitConvert[pos, "Degrees"];
GST = UnitConvert[SiderealTime[pos0, date0], "Degrees"];
r = PlanetData["Earth",
EntityProperty["Planet", "DistanceFromSun", {"Date" -> date0}]];
long = Quantity[360, "Degrees"] - (GST - pos[[1]]);
lat = pos[[2]];
<|
x -> r Cos[lat] Cos[long],
y -> r Cos[lat] Sin[long],
z -> r Sin[lat]
|>
]
]
Moon position
MoonPositionXYZ[date_] := With[
{date0 = TimeZoneConvert[date, 0], pos0 = GeoPosition[{0, 0}]},
Module[{pos, GST, r, lat, long},
pos = PlanetaryMoonData["Moon",
{EntityProperty["PlanetaryMoon", "RightAscension", {"Date" -> date0}],
EntityProperty["PlanetaryMoon", "Declination", {"Date" -> date0}]}];
pos = UnitConvert[pos, "Degrees"];
GST = UnitConvert[SiderealTime[pos0, date0], "Degrees"];
r = PlanetaryMoonData["Moon",
EntityProperty["PlanetaryMoon", "DistanceFromSun", {"Date" -> date0}]];
long = Quantity[360, "Degrees"] - (GST - pos[[1]]);
lat = pos[[2]];
<|
x -> r Cos[lat] Cos[long],
y -> r Cos[lat] Sin[long],
z -> r Sin[lat]
|>
]
]
Planet positions
PlanetPositionXYZ[planet_, date_] := With[
{date0 = TimeZoneConvert[date, 0], pos0 = GeoPosition[{0, 0}]},
Module[{pos, GST, r, lat, long},
pos = PlanetData[planet,
{EntityProperty["Planet", "RightAscension", {"Date" -> date0}],
EntityProperty["Planet", "Declination", {"Date" -> date0}]}];
pos = UnitConvert[pos, "Degrees"];
GST = UnitConvert[SiderealTime[pos0, date0], "Degrees"];
r = PlanetData[planet,
EntityProperty["Planet", "DistanceFromEarth", {"Date" -> date0}]];
long = Quantity[360, "Degrees"] - (GST - pos[[1]]);
lat = pos[[2]];
<|
x -> r Cos[lat] Cos[long],
y -> r Cos[lat] Sin[long],
z -> r Sin[lat]
|>
]
]
For moons, just make the substitution PlanetData -> PlanetaryMoonData and "Planet" -> "PlanetaryMoon".
Date and Time
Run these using DateObject to specify the time and date:
MoonPositionXYZ[DateObject[{2015, 4, 17, 10, 0, 0}, TimeZone -> -3]]
SunPositionXYZ[DateObject[{2015, 4, 17, 10, 0, 0}, TimeZone -> 10]]
PlanetPositionXYZ["Venus",
DateObject[{2015, 4, 17, 10, 0, 0}, TimeZone -> 5]]
(* -> <|x -> 0.857983 au, y -> -0.512016 au, z -> 0.0685721 au|>
<|x -> -0.9875 au, y -> -0.00628945 au, z -> 0.179209 au|>
<|x -> -0.824832 au, y -> 0.583867 au, z -> 0.431657 au|> *)
Mathematica should convert all of the DateObjects into some kind of universal date (probably Julian date) internally when they are passed to the various astronomical position functions, so as long as the date and time are uniquely specified, it shouldn't matter what timezone is used (note that Mathematica may use your timezone as default). That said, it is probably safest to specify the time and date fully (as above), rather than just using DateObject[TimeZone -> x], as this seems to give an unexpected date sometimes (see this question).
Other objects
For objects such as stars for which Mathematica gives non-precessing RA and Dec, we must first convert it to a precessing RA. Mathematica already does this, but the functionality is buried in the innards of the astronomical Data functions. I have not yet come up with a good way of doing the conversion.
