# skyfield

## How to compute & plot sun path diagram using skyfield in Python

In this example, we’ll show how to generate a sun path diagram for any given day for a preselected location. By using the skyfield library, we can obtain extremely accurate predictions of the sun’s position in the sky, and in contrast to many other libraries, we can use the same method to predict the position of other celestial bodies.

```#!/usr/bin/env python3
from skyfield import api
from skyfield import almanac
from datetime import datetime
from datetime import timedelta
from matplotlib import pyplot as plt
import pytz
import dateutil.parser
from collections import namedtuple
from UliEngineering.Utils.Date import *
import matplotlib.ticker as mtick

# Matplotlib formatter
def format_degrees(value, pos=None):
return f'{value:.0f} °'

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

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

# Data types
AltAz = namedtuple("AltAz", ["altitude", "azimuth", "distance"])
TimeAltAz = namedtuple("TimeAltAz", ["time", "altitude", "azimuth", "distance"])

def sun_altaz(planet_at_location, ts):
# Compute the sun position as seen from the observer at <location>
sun_pos = planet_at_location.at(ts).observe(sun).apparent()
# Compute apparent altitude & azimuth for the sun's position
altitude, azimuth, distance = sun_pos.altaz()
return AltAz(altitude, azimuth, distance)

def sun_altaz_for_day(location, year, month, day, tz):
earth_at_location = (earth + location)

minutes = list(yield_minutes_on_day(year=year, month=month, day=day, tz=tz))
skyfield_minutes = ts.from_datetimes(minutes)

minutely_altaz = [sun_altaz(earth_at_location, ts) for ts in skyfield_minutes]
# Extract components for plotting
minutely_alt = [altaz.altitude.degrees for altaz in minutely_altaz]
minutely_az = [altaz.azimuth.degrees for altaz in minutely_altaz]
minutely_distance = [altaz.distance for altaz in minutely_altaz]

return TimeAltAz(minutes, minutely_alt, minutely_az, minutely_distance)

# Random location near Munich
location = api.Topos('48.324777 N', '11.405610 E', elevation_m=519)
# Compute Alt/Az for two different days
jun = sun_altaz_for_day(location, 2022, 6, 15, pytz.timezone("Europe/Berlin"))
dec = sun_altaz_for_day(location, 2022, 12, 15, pytz.timezone("Europe/Berlin"))

# Plot!
plt.gca().yaxis.set_major_formatter(mtick.FuncFormatter(format_degrees))
plt.gca().xaxis.set_major_formatter(mtick.FuncFormatter(format_degrees))

plt.plot(jun.azimuth, jun.altitude, label="15th of June")
plt.plot(dec.azimuth, dec.altitude, label="15th of December")
plt.ylim([0, None]) # Do not show sun angles below the horizon
plt.ylabel("Sun altitude")
plt.xlabel("Sun azimuth")
plt.title("Apparent sun position at our office")
plt.grid()
plt.gca().legend()
plt.gcf().set_size_inches(10,5)

plt.savefig("Sun path.svg")```

Posted by Uli Köhler in Physics, Python, skyfield

## How to convert skyfield Time into datetime at specific timezone

When you have a skyfield `Time` object like

```t = ts.now()
# Example: <Time tt=2459750.027604357>
```

you can convert it to a Python `datetime` in a specific timezone (`Europe/Berlin` in this example) using `.astimezone()` and the `pytz` library:

```t.astimezone(pytz.timezone("Europe/Berlin"))
# Example: datetime.datetime(2022, 6, 19, 14, 38, 35, 832445, tzinfo=<DstTzInfo 'Europe/Berlin' CEST+2:00:00 DST>)```

### Complete example

```from skyfield import api
from datetime import datetime
import pytz

ts = api.load.timescale()
t = ts.now()
dt = t.astimezone(pytz.timezone("Europe/Berlin"))
print(dt) # e.g. 2022-06-19 14:42:47.406786+02:00
```

Posted by Uli Köhler in Python, skyfield

## 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
Posted by Uli Köhler in Physics, Python, skyfield

## How to fix Python skyfield FileNotFoundError: [Errno 2] No such file or directory: ‘de421.bsp’

### Problem:

When trying to use the Python `skyfield` library, you see an exception like

```Input In [2], in <cell line: 11>()
8 from calendar import monthrange
10 ts = api.load.timescale()
---> 11 ephem = api.load_file('de413.bsp')

File /usr/local/lib/python3.10/dist-packages/skyfield/iokit.py:412, in load_file(path)
410 base, ext = os.path.splitext(path)
411 if ext == '.bsp':
--> 412     return SpiceKernel(path)
413 raise ValueError('unrecognized file extension: {}'.format(path))

File /usr/local/lib/python3.10/dist-packages/skyfield/jpllib.py:71, in SpiceKernel.__init__(self, path)
69 self.path = path
70 self.filename = os.path.basename(path)
---> 71 self.spk = SPK.open(path)
72 self.segments = [SPICESegment(self, s) for s in self.spk.segments]
73 self.codes = set(s.center for s in self.segments).union(
74                  s.target for s in self.segments)

File /usr/local/lib/python3.10/dist-packages/jplephem/spk.py:49, in SPK.open(cls, path)
46 @classmethod
47 def open(cls, path):
48     """Open the file at `path` and return an SPK instance."""
---> 49     return cls(DAF(open(path, 'rb')))

FileNotFoundError: [Errno 2] No such file or directory: 'de421.bsp'```

### Solution:

Take a look at the `api.load(...)` line in your code:

`ephem = api.load_file('de421.bsp')`

It tries to load the data from the file `de421.bsp` in the current directory. This file contains positional data of objects in the sky and you need to manually download that file.

You can download the file from NASA. Just take care to either place it into the right directory or modifying the path in the `api.load()` call to point to the file.

URL for downloading the file:

`https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/de421.bsp`

My preferred way to download it is using `wget`:

`wget https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/de421.bsp`

This command will place the file into the current directory.

Posted by Uli Köhler in Python, skyfield

## How to compute sunrise & sunset in Python using skyfield

The following code will compute the sunrise & sunset at a specific location & elevation using the `skyfield` library. Note that you need to download `de413.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('de413.bsp')

def compute_sunrise_sunset(location, year=2019, month=1, day=1):
t0 = ts.utc(year, month, day, 0)
# t1 = t0 plus one day
t1 = ts.utc(t0.utc_datetime() + timedelta(days=1))
t, y = almanac.find_discrete(t0, t1, almanac.sunrise_sunset(ephem, location))
sunrise = None
for time, is_sunrise in zip(t, y):
if is_sunrise:
sunrise = dateutil.parser.parse(time.utc_iso())
else:
sunset = dateutil.parser.parse(time.utc_iso())
return sunrise, sunset

# Compute sunrise & sunset for random location near Munich
location = api.Topos('48.324777 N', '11.405610 E', elevation_m=519)
now = datetime.now()
sunrise, sunset = compute_sunrise_sunset(location, now.year, now.month, now.day)

# Print result (example)
print(f'Sunrise today: {sunrise}')
print(f'Sunset today: {sunset}')```

#### Definition of sunrise & sunset in this context

According to the `skyfield` documentation:

`Skyfield uses the same definition as the United States Naval Observatory: the Sun is up when its center is 0.8333 degrees below the horizon, which accounts for both its apparent radius of around 16 arcminutes and also for the 34 arcminutes by which atmospheric refraction on average lifts the image of the Sun.`

#### Other caveats

• Note that obstructions like mountains are not taken into account for this model
• Note that the resulting timestamps are UTC, if you want local time, you’ll have to convert them appropriately

#### Example output:

```Sunrise today: 2022-06-19 03:12:56+00:00
Sunset today: 2022-06-19 19:18:38+00:00```

Posted by Uli Köhler in Physics, Python, skyfield

## How to fix Python skyfield FileNotFoundError: [Errno 2] No such file or directory: ‘de413.bsp’

### Problem:

When trying to use the Python `skyfield` library, you see an exception like

```Input In [2], in <cell line: 11>()
8 from calendar import monthrange
10 ts = api.load.timescale()
---> 11 ephem = api.load_file('de413.bsp')

File /usr/local/lib/python3.10/dist-packages/skyfield/iokit.py:412, in load_file(path)
410 base, ext = os.path.splitext(path)
411 if ext == '.bsp':
--> 412     return SpiceKernel(path)
413 raise ValueError('unrecognized file extension: {}'.format(path))

File /usr/local/lib/python3.10/dist-packages/skyfield/jpllib.py:71, in SpiceKernel.__init__(self, path)
69 self.path = path
70 self.filename = os.path.basename(path)
---> 71 self.spk = SPK.open(path)
72 self.segments = [SPICESegment(self, s) for s in self.spk.segments]
73 self.codes = set(s.center for s in self.segments).union(
74                  s.target for s in self.segments)

File /usr/local/lib/python3.10/dist-packages/jplephem/spk.py:49, in SPK.open(cls, path)
46 @classmethod
47 def open(cls, path):
48     """Open the file at `path` and return an SPK instance."""
---> 49     return cls(DAF(open(path, 'rb')))

FileNotFoundError: [Errno 2] No such file or directory: 'de413.bsp'```

### Solution:

Take a look at the `api.load(...)` line in your code:

`ephem = api.load_file('de413.bsp')`

It tries to load the data from the file `de413.bsp` in the current directory. This file contains positional data of objects in the sky and you need to manually download that file.

You can download the file from NASA. Just take care to either place it into the right directory or modifying the path in the `api.load()` call to point to the file.

URL for downloading the file:

`https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/de413.bsp`

My preferred way to download it is using `wget`:

`wget https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/de413.bsp`

This command will place the file into the current directory.

Posted by Uli Köhler in Python, skyfield