<p>I'm trying to use astropy to plot a simulation of the solar eclipse that we had in Exmouth 20 April 2023 11:31 Local time. I'm a programmer not a physist! So after much messing around I tried this:</p>
<pre><code>import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from matplotlib.animation import FuncAnimation
from astropy.coordinates import EarthLocation, AltAz, get_body,solar_system_ephemeris
from astropy.time import Time
from astropy import units as u
import os
from matplotlib.transforms import Bbox
#SERRURIER ISLAND 21.6087° S, 114.6802° E
# Define the location of Exmouth
with solar_system_ephemeris.set('jpl'):
    location = EarthLocation.from_geodetic(114.6802*u.deg,-21.6087*u.deg,0*u.m)

    # Define the time range for the eclipse
    start_time = Time('2023-04-20 02:00:00')  # Start time (UTC)
    end_time = Time('2023-04-20 06:00:00')    # End time (UTC)
    time_resolution = 60 * u.second             # Time resolution (1 second)
    times = start_time + np.arange(0, (end_time - start_time).to(u.second).value + 1, time_resolution.value) * u.second

    # Convert times to local time in Exmouth (UTC+8)
    local_times = times + 8*u.hour

    # Get the positions of the sun and moon
    alt_az = AltAz(obstime=times, location=location)

    sun_altaz = get_body('sun',times,ephemeris='de440s').transform_to(alt_az)
    moon_altaz = get_body('moon',times,ephemeris='de440s').transform_to(alt_az)
    sun_moon_sep = moon_altaz.separation(sun_altaz)
    i0 = sun_moon_sep.argmin() - 128
    i1 = sun_moon_sep.argmin() + 128
    # Calculate the angular diameters
    sun_diameter = 1_391_000 * u.km
    moon_diameter = 3_474 * u.km
    distance_to_sun = 149_600_000 * u.km
    distance_to_moon = 384_400 * u.km

    sun_angular_radius = ((sun_diameter / sun_altaz.distance)*u.rad).to(u.deg).value
    moon_angular_radius = ((moon_diameter / moon_altaz.distance)*u.rad).to(u.deg).value
    # Calculate the angular sizes of the sun and moon
    #sun_angular_radius = np.arctan2(696340*u.km, get_body('sun',times).distance).to(u.deg)
    #moon_angular_radius = np.arctan2(1737.4*u.km, get_body('moon',times).distance).to(u.deg)

# Create a directory to save the frames
output_dir = '/home/mor582/eclipse_frames'
os.makedirs(output_dir, exist_ok=True)
# Set up the plot
fig, ax = plt.subplots(figsize=(19.2,10.8), dpi=100)
plt.gca().set_position([0, 0, 1, 1])


ax.set_xlim(-8, 8)
#ax.set_ylim(-8*0.5625, 8*0.5625)
ax.set_aspect('equal')
ax.set_xlabel('Degrees')
ax.set_ylabel('Degrees')
ax.set_axis_off()
xlim = ax.get_xlim()
ylim = ax.get_ylim()
# Plot the sun and moon

moon_circle = mpl.patches.Ellipse((0,0), width=moon_angular_radius[0], height=moon_angular_radius[0],
                              color='#666666', zorder=10)
sun_circle = mpl.patches.Ellipse((0,0), width=sun_angular_radius[0], height=sun_angular_radius[0],
                             color='orange', zorder=1)
# sun_circle = plt.Circle((0.0, 0.0), sun_angular_radius.value, color='yellow', ec='black', lw=2)
# # Set the new position for the time text
#moon_circle = plt.Circle((0.0, 0.0), moon_angular_radius.value, color='gray', ec='black', lw=2)
ax.add_artist(sun_circle)
ax.add_artist(moon_circle)

# Add text annotation for the local time with a semi-transparent background
time_text = ax.text(0, 0, '2023-04-20 11:00:00', fontsize=12, ha='center', color='black',
                    bbox=dict(facecolor='white', alpha=0.6, edgecolor='none'),fontweight='bold')
time_text_bbox = time_text.get_window_extent(renderer=fig.canvas.get_renderer(),dpi=100)
bbox_text = Bbox(ax.transData.inverted().transform(time_text_bbox))
time_text_x = xlim[1] - 0.1 -bbox_text.width/2 
time_text_y = ylim[1] -0.01 -bbox_text.height 

# Set the new position for the time text
time_text.set_position((time_text_x, time_text_y))

sun_bbox = sun_circle.get_window_extent(renderer=fig.canvas.get_renderer())
bbox_sun = Bbox(ax.transData.inverted().transform(sun_bbox))
sun_x = xlim[1] - bbox_sun.width/2 
sun_y = time_text_y -ax.margins()[1] -bbox_sun.height
#sun_circle.center =(sun_x,sun_y)


# Set the background color to be transparent
fig.patch.set_alpha(0)
ax.patch.set_alpha(0)
ax.axis('off')

def wrap_north(bearing):
    if bearing&gt;180:
        bearing = bearing -360
    return bearing

# Function to update the position of the moon and the sun's transparency in the plot
def update(frame):
    # Update the moon's position
    moon_circle.center = (wrap_north(moon_altaz.az[frame].deg),moon_altaz.alt[frame].deg)
    sun_circle.center = (wrap_north(sun_altaz.az[frame].deg),sun_altaz.alt[frame].deg)
    az_lim = (wrap_north(sun_altaz.az[frame].deg) - 2,
              wrap_north(sun_altaz.az[frame].deg) + 2)
    #alt_lim = (sun_aa.alt[i].to(u.degree).value - 2,
    #           sun_aa.alt[i].to(u.degree).value + 2)
    alt_lim = (sun_altaz.alt[frame].deg - 2,
               sun_altaz.alt[frame].deg + 2)
    ax.set_xlim(az_lim)
    ax.set_ylim(alt_lim)
    time_text_bbox = time_text.get_window_extent(renderer=fig.canvas.get_renderer(),dpi=100)
    bbox_text = Bbox(ax.transData.inverted().transform(time_text_bbox))
    time_text_x = az_lim[1]-bbox_text.width/2 
    time_text_y = alt_lim[1] -bbox_text.height *2
    #moon_circle.height = moon_angular_radius * np.cos(moon_altaz.alt[frame])
    #sun_circle.height = sun_angular_radius * np.cos(sun_altaz.alt[frame])
    # Calculate the alpha value for the sun's transparency based on its altitude
    #alpha = max(0, min(1, (sun_altitude - sun_altaz.alt.min().deg) / (sun_altaz.alt.max().deg - sun_altaz.alt.min().deg)))
    #sun_circle.set_alpha(alpha)
    
    # Update the local time text

    local_time_str = local_times[frame].iso[:-4]  # Format time string without seconds
    time_text.set_text(local_time_str)
    time_text.set_position([time_text_x,time_text_y])

    # Position the time text in the top-right corner
    
    # Save the frame as a PNG file with a transparent background
    file_time=local_times[frame].iso.replace('-', '').replace(':', '').replace(' ', 'T').split('.')[0]
    hour_dir =f&quot;{output_dir}/{file_time[:-4]}&quot;
    # if not os.path.exists(hour_dir):
    #     os.mkdir(hour_dir)
    #/{file_time[:-4]}
    plt.savefig(f&quot;{output_dir}/frame_{file_time}_{frame:04d}.png&quot;, transparent=True)
    return sun_circle, moon_circle, time_text

# Create animation
ani = FuncAnimation(fig, update, frames=len(times), blit=True, interval=1000)

# Save each frame as an image
for i in range(len(times)):
    update(i)

plt.close()
</code></pre>
<p>Here is a link to what I'm trying to recreate</p>
<p><a href="https://www.timeanddate.com/eclipse/in/australia/exmouth?iso=20230420" rel="nofollow noreferrer">https://www.timeanddate.com/eclipse/in/australia/exmouth?iso=20230420</a></p>
<p>There seems to be a difference in the Sun's position between Astropy and this rendering.</p>
<p><a href="https://i.sstatic.net/zf7bpT5n.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/zf7bpT5n.png" alt="enter image description here" /></a></p>
<p>If I use sunpy and the sun.eclipse_amount(observer, moon_radius='minimum') it gives the correct times for the start and end of the eclipse.</p>

CC BY-SA 4.0
2024-10-28T10:19:51.873Z
Still looking, but one issue might be that you've hard-coded the distances to the sun and moon, you should use the values stored in the SkyCoord instance.
Roy Smart
CC BY-SA 4.0
2024-10-28T11:57:36.537Z
Thanks yep I don't understand using the sunpy amount function for the angular separation of the moon and the sun  about 0.5 at 13:04 local but they are using a different coordinate system
Nick Mortimer
CC BY-SA 4.0
2024-10-28T12:07:13.863Z
Here's some geocentric data from JPL Horizons. Sun https://ssd.jpl.nasa.gov/api/horizons.api?format=text&COMMAND=10&CENTER=%40399&START_TIME=%272923-Apr-20+11%3A30+UT%2B8%27&STOP_TIME=%272923-Apr-20+11%3A32%27&STEP_SIZE=1m&QUANTITIES=%2720%2C13%27&RANGE_UNITS=KM&OBJ_DATA=NO Moon https://ssd.jpl.nasa.gov/api/horizons.api?format=text&COMMAND=301&CENTER=%40399&START_TIME=%272923-Apr-20+11%3A30+UT%2B8%27&STOP_TIME=%272923-Apr-20+11%3A32%27&STEP_SIZE=1m&QUANTITIES=%2720%2C13%27&RANGE_UNITS=KM&OBJ_DATA=NO
PM 2Ring
CC BY-SA 4.0
2024-10-28T16:20:53.310Z
Usually Besselian Elements are used to compute the paths and local circumstances.  Here is an open source implementation, references are in the eclipse.js file.  https://celestialprogramming.com/apps/SolarEclipseViewer/viewer.html
Greg Miller
CC BY-SA 4.0
2024-10-29T02:29:18.507Z
Thanks Greg, I'm trying to create a graphic of the moon passing over the sun as seen at Exmouth Western Australia rather than the track of the eclipse. The function amount in sunpy correctly calculates the angular separation but I don't know how to plot it to give the desired result
Nick Mortimer
<astropy><eclipse>
Plotting the progress of an eclipse
Nick Mortimer
# Plotting the progress of an eclipse

(Asked by: username_0 on 2024-10-28)

I'm trying to use astropy to plot a simulation of the solar eclipse that we had in Exmouth 20 April 2023 11:31 Local time. I'm a programmer not a physist! So after much messing around I tried this:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from matplotlib.animation import FuncAnimation
from astropy.coordinates import EarthLocation, AltAz, get_body,solar_system_ephemeris
from astropy.time import Time
from astropy import units as u
import os
from matplotlib.transforms import Bbox
#SERRURIER ISLAND 21.6087° S, 114.6802° E
# Define the location of Exmouth
with solar_system_ephemeris.set('jpl'):
    location = EarthLocation.from_geodetic(114.6802*u.deg,-21.6087*u.deg,0*u.m)

    # Define the time range for the eclipse
    start_time = Time('2023-04-20 02:00:00')  # Start time (UTC)
    end_time = Time('2023-04-20 06:00:00')    # End time (UTC)
    time_resolution = 60 * u.second             # Time resolution (1 second)
    times = start_time + np.arange(0, (end_time - start_time).to(u.second).value + 1, time_resolution.value) * u.second

    # Convert times to local time in Exmouth (UTC+8)
    local_times = times + 8*u.hour

    # Get the positions of the sun and moon
    alt_az = AltAz(obstime=times, location=location)

    sun_altaz = get_body('sun',times,ephemeris='de440s').transform_to(alt_az)
    moon_altaz = get_body('moon',times,ephemeris='de440s').transform_to(alt_az)
    sun_moon_sep = moon_altaz.separation(sun_altaz)
    i0 = sun_moon_sep.argmin() - 128
    i1 = sun_moon_sep.argmin() + 128
    # Calculate the angular diameters
    sun_diameter = 1_391_000 * u.km
    moon_diameter = 3_474 * u.km
    distance_to_sun = 149_600_000 * u.km
    distance_to_moon = 384_400 * u.km

    sun_angular_radius = ((sun_diameter / sun_altaz.distance)*u.rad).to(u.deg).value
    moon_angular_radius = ((moon_diameter / moon_altaz.distance)*u.rad).to(u.deg).value
    # Calculate the angular sizes of the sun and moon
    #sun_angular_radius = np.arctan2(696340*u.km, get_body('sun',times).distance).to(u.deg)
    #moon_angular_radius = np.arctan2(1737.4*u.km, get_body('moon',times).distance).to(u.deg)

# Create a directory to save the frames
output_dir = '/home/mor582/eclipse_frames'
os.makedirs(output_dir, exist_ok=True)
# Set up the plot
fig, ax = plt.subplots(figsize=(19.2,10.8), dpi=100)
plt.gca().set_position([0, 0, 1, 1])


ax.set_xlim(-8, 8)
#ax.set_ylim(-8*0.5625, 8*0.5625)
ax.set_aspect('equal')
ax.set_xlabel('Degrees')
ax.set_ylabel('Degrees')
ax.set_axis_off()
xlim = ax.get_xlim()
ylim = ax.get_ylim()
# Plot the sun and moon

moon_circle = mpl.patches.Ellipse((0,0), width=moon_angular_radius[0], height=moon_angular_radius[0],
                              color='#666666', zorder=10)
sun_circle = mpl.patches.Ellipse((0,0), width=sun_angular_radius[0], height=sun_angular_radius[0],
                             color='orange', zorder=1)
# sun_circle = plt.Circle((0.0, 0.0), sun_angular_radius.value, color='yellow', ec='black', lw=2)
# # Set the new position for the time text
#moon_circle = plt.Circle((0.0, 0.0), moon_angular_radius.value, color='gray', ec='black', lw=2)
ax.add_artist(sun_circle)
ax.add_artist(moon_circle)

# Add text annotation for the local time with a semi-transparent background
time_text = ax.text(0, 0, '2023-04-20 11:00:00', fontsize=12, ha='center', color='black',
                    bbox=dict(facecolor='white', alpha=0.6, edgecolor='none'),fontweight='bold')
time_text_bbox = time_text.get_window_extent(renderer=fig.canvas.get_renderer(),dpi=100)
bbox_text = Bbox(ax.transData.inverted().transform(time_text_bbox))
time_text_x = xlim[1] - 0.1 -bbox_text.width/2 
time_text_y = ylim[1] -0.01 -bbox_text.height 

# Set the new position for the time text
time_text.set_position((time_text_x, time_text_y))

sun_bbox = sun_circle.get_window_extent(renderer=fig.canvas.get_renderer())
bbox_sun = Bbox(ax.transData.inverted().transform(sun_bbox))
sun_x = xlim[1] - bbox_sun.width/2 
sun_y = time_text_y -ax.margins()[1] -bbox_sun.height
#sun_circle.center =(sun_x,sun_y)


# Set the background color to be transparent
fig.patch.set_alpha(0)
ax.patch.set_alpha(0)
ax.axis('off')

def wrap_north(bearing):
    if bearing>180:
        bearing = bearing -360
    return bearing

# Function to update the position of the moon and the sun's transparency in the plot
def update(frame):
    # Update the moon's position
    moon_circle.center = (wrap_north(moon_altaz.az[frame].deg),moon_altaz.alt[frame].deg)
    sun_circle.center = (wrap_north(sun_altaz.az[frame].deg),sun_altaz.alt[frame].deg)
    az_lim = (wrap_north(sun_altaz.az[frame].deg) - 2,
              wrap_north(sun_altaz.az[frame].deg) + 2)
    #alt_lim = (sun_aa.alt[i].to(u.degree).value - 2,
    #           sun_aa.alt[i].to(u.degree).value + 2)
    alt_lim = (sun_altaz.alt[frame].deg - 2,
               sun_altaz.alt[frame].deg + 2)
    ax.set_xlim(az_lim)
    ax.set_ylim(alt_lim)
    time_text_bbox = time_text.get_window_extent(renderer=fig.canvas.get_renderer(),dpi=100)
    bbox_text = Bbox(ax.transData.inverted().transform(time_text_bbox))
    time_text_x = az_lim[1]-bbox_text.width/2 
    time_text_y = alt_lim[1] -bbox_text.height *2
    #moon_circle.height = moon_angular_radius * np.cos(moon_altaz.alt[frame])
    #sun_circle.height = sun_angular_radius * np.cos(sun_altaz.alt[frame])
    # Calculate the alpha value for the sun's transparency based on its altitude
    #alpha = max(0, min(1, (sun_altitude - sun_altaz.alt.min().deg) / (sun_altaz.alt.max().deg - sun_altaz.alt.min().deg)))
    #sun_circle.set_alpha(alpha)
    
    # Update the local time text

    local_time_str = local_times[frame].iso[:-4]  # Format time string without seconds
    time_text.set_text(local_time_str)
    time_text.set_position([time_text_x,time_text_y])

    # Position the time text in the top-right corner
    
    # Save the frame as a PNG file with a transparent background
    file_time=local_times[frame].iso.replace('-', '').replace(':', '').replace(' ', 'T').split('.')[0]
    hour_dir =f"{output_dir}/{file_time[:-4]}"
    # if not os.path.exists(hour_dir):
    #     os.mkdir(hour_dir)
    #/{file_time[:-4]}
    plt.savefig(f"{output_dir}/frame_{file_time}_{frame:04d}.png", transparent=True)
    return sun_circle, moon_circle, time_text

# Create animation
ani = FuncAnimation(fig, update, frames=len(times), blit=True, interval=1000)

# Save each frame as an image
for i in range(len(times)):
    update(i)

plt.close()

Here is a link to what I'm trying to recreate
https://www.timeanddate.com/eclipse/in/australia/exmouth?iso=20230420
There seems to be a difference in the Sun's position between Astropy and this rendering.

If I use sunpy and the sun.eclipse_amount(observer, moon_radius='minimum') it gives the correct times for the start and end of the eclipse.

**Tags:** astropy, eclipse
> **Comments:**
> **username_1**: Still looking, but one issue might be that you've hard-coded the distances to the sun and moon, you should use the values stored in the SkyCoord instance.  
> **username_0**: Thanks yep I don't understand using the sunpy amount function for the angular separation of the moon and the sun  about 0.5 at 13:04 local but they are using a different coordinate system  
> **username_2**: Here's some geocentric data from JPL Horizons. Sun https://ssd.jpl.nasa.gov/api/horizons.api?format=text&COMMAND=10&CENTER=%40399&START_TIME=%272923-Apr-20+11%3A30+UT%2B8%27&STOP_TIME=%272923-Apr-20+11%3A32%27&STEP_SIZE=1m&QUANTITIES=%2720%2C13%27&RANGE_UNITS=KM&OBJ_DATA=NO Moon https://ssd.jpl.nasa.gov/api/horizons.api?format=text&COMMAND=301&CENTER=%40399&START_TIME=%272923-Apr-20+11%3A30+UT%2B8%27&STOP_TIME=%272923-Apr-20+11%3A32%27&STEP_SIZE=1m&QUANTITIES=%2720%2C13%27&RANGE_UNITS=KM&OBJ_DATA=NO  
> **username_3**: Usually Besselian Elements are used to compute the paths and local circumstances.  Here is an open source implementation, references are in the eclipse.js file.  https://celestialprogramming.com/apps/SolarEclipseViewer/viewer.html  
> **username_0**: Thanks Greg, I'm trying to create a graphic of the moon passing over the sun as seen at Exmouth Western Australia rather than the track of the eclipse. The function amount in sunpy correctly calculates the angular separation but I don't know how to plot it to give the desired result  

---

## 0 answers

