Absolutely! Skyfield is more than capable of calculating both the sub-Earth point (the spot on the Moon directly facing Earth) and sub-Solar point (the spot directly facing the Sun) for your lunar crater height project. Here’s a straightforward, actionable implementation:

First, make sure you have Skyfield installed (if not, run pip install skyfield), then follow these steps:

Step 1: Set Up Dependencies & Celestial Bodies #

Start by importing the necessary tools and loading a high-precision ephemeris (Skyfield will automatically download the required data if you don’t have it locally):

from skyfield.api import load, utc

# Initialize timescale and ephemeris
ts = load.timescale()
eph = load('de421.bsp')  # DE421 is a widely used ephemeris for lunar/solar system calculations

# Define key celestial bodies
earth = eph['earth']
moon = eph['moon']
sun = eph['sun']

Step 2: Specify Your Observation Time #

You’ll need the UTC timestamp when you captured your lunar image (convert your local time to UTC for accuracy—this is critical for precise calculations):

# Replace with your actual image capture time (UTC)
observation_time = ts.utc(2024, 5, 20, 18, 45, 0)  # Example: 2024-05-20 18:45:00 UTC

Step 3: Calculate the Sub-Earth Point #

The sub-Earth point is the location on the Moon where Earth is directly overhead. To find this, we calculate Earth’s position relative to the Moon, then convert it to the Moon’s fixed coordinate system (iau_moon):

# Get Earth's position relative to the Moon at the observation time
earth_relative_to_moon = moon.at(observation_time).observe(earth)

# Convert to Moon's fixed (IAU) coordinate system to get lat/lon
sub_earth_lat, sub_earth_lon, _ = earth_relative_to_moon.frame_xyz('iau_moon').latlon()

# Print results in degrees
print(f"Sub-Earth Point:")
print(f"  Latitude: {sub_earth_lat.degrees:.2f}°")
print(f"  Longitude: {sub_earth_lon.degrees:.2f}°")

Step 4: Calculate the Sub-Solar Point #

Similarly, the sub-Solar point is where the Sun is directly overhead on the Moon. We use the same approach but with the Sun’s position relative to the Moon:

# Get Sun's position relative to the Moon at the observation time
sun_relative_to_moon = moon.at(observation_time).observe(sun)

# Convert to Moon's fixed coordinate system
sub_solar_lat, sub_solar_lon, _ = sun_relative_to_moon.frame_xyz('iau_moon').latlon()

# Print results in degrees
print(f"\nSub-Solar Point:")
print(f"  Latitude: {sub_solar_lat.degrees:.2f}°")
print(f"  Longitude: {sub_solar_lon.degrees:.2f}°")

Key Notes for Your Project #

• The iau_moon coordinate system uses the Moon’s official fixed reference frame (aligned with its rotational poles), which is exactly what you need for consistent crater position calculations.
• Double-check your UTC timestamp—even a 10-minute error can shift the sub-point positions enough to affect your crater height measurements.
• If you need higher precision, you can use a more recent ephemeris like de440.bsp (just replace de421.bsp in the load() call).

This implementation should give you the exact coordinates required by your lunar tutorial, and it’s fully compatible with Skyfield’s modern, Pythonic API (no need to rely on PyEphem anymore).

内容的提问来源于stack exchange，提问作者Dimitrios Theodorakis