How to compute position of sun in the sky in Python using skyfield

The following code computes the position of the sun (in azimuth/altitude coordinates) in the sky using the skyfield library. Note that you need to download de421.bsp .

from skyfield import api
from skyfield import almanac
from datetime import datetime
from datetime import timedelta
import dateutil.parser
from calendar import monthrange

ts = api.load.timescale()
ephem = api.load_file('de421.bsp')

sun = ephem["Sun"]
earth = ephem["Earth"]

# Compute sunrise & sunset for random location near Munich
location = api.Topos('48.324777 N', '11.405610 E', elevation_m=519)
# Compute the sun position as seen from the observer at <location>
sun_pos = (earth + location).at(ts.now()).observe(sun).apparent()
# Compute apparent altitude & azimuth for the sun's position
altitude, azimuth, distance = sun_pos.altaz()

# Print results (example)
print(f"Altitude: {altitude.degrees:.4f} °")
print(f"Azimuth: {azimuth.degrees:.4f} °")

Example output

Altitude: -3.3121 °
Azimuth: 48.4141 °

Factors influencing the accuracy of the calculation

This way of calculating the position takes into account:

The slight shift in position caused by light speed
The very very slight shift in position caused by earth’s gravity

But it does not take into account:

Atmospheric distortions shifting the sun’s position
The extent of the sun’s disk causing the sun to emanate not from a point but apparently from an area