I Made a Sky Map in Python. Here’s How.

This tutorial shows you how to create a sky map / chart in Python using skyfield library.

Photo by Andy Holmes on Unsplash

“For my part I know nothing with any certainty, but the sight of the stars makes me dream.” — Vincent Van Gogh

There has always been something so fascinating about stars: the way they shine, align into constellations, hold so much unknown… Among many things, humans have been interested in mapping out stars since the dawn of time. Star charts have been used primarily for navigation but are also useful for various scientific and exploratory purposes. Celestial maps have also made their way into our homes in form of custom art prints or digital downloads (as evident by 15K+ results on Etsy).

There are numerous websites that can generate star charts for free or for a fee but I really wanted to learn how to do it myself so I can control many aspects of it. And here, I’d like to share what I learned with you. Let’s get started.

Prerequisites

First, you need a version of Python installed on your machine. I am working out of Python 3.8.8 on my Mac Big .

Then, we need the following libraries. I have listed them with the specific version I used in the demo.

• Numpy 1.19.2 — used for working with arrays of star data
• Matplotlib 3.3.4 — used for making the final star map image
• Skyfield 1.45 — used for getting star data and calculating star positions
• Datetime — used for converting string dates into date-time objects

The following libraries are optional but useful if like me, you are also getting old, lazy, and cannot be bothered:

• Geopy 2.2.0 — converts our location into latitude and longitude
• Tzwhere — grabs the timezone of our latitude and longitude
• Pytz 2021.3 — used for converting date-time from local timezone into UTC timezone (which is needed for Skyfield to build the sky map properly)

Now that we have installed and upgraded everything we need, let’s import our libraries, objects and methods and get started.

# datetime libraries
from datetime import datetime
from geopy import Nominatim
from tzwhere import tzwhere
from pytz import timezone, utc
# matplotlib to help display our star map
import matplotlib.pyplot as plt
# skyfield for star data 
from skyfield.api import Star, load, wgs84
from skyfield.data import hipparcos
from skyfield.projections import build_stereographic_projection

Loading Earth and Star Data

First, we need to use the skyfield library in order to load our star data.

# de421 shows position of earth and sun in space
eph = load('de421.bsp')
# hipparcos dataset contains star location data
with load.open(hipparcos.URL) as f:
    stars = hipparcos.load_dataframe(f)

Setting Up the Observation Location

In the first programming step, we need to define the date and time of when we want to capture our sky. We also want to define the location as a skyfield object and use it to make star location calculation later on.

As an example for this demo, I want to get a star map for Times Square in NYC on New Year’s Eve. I am going to define these parameters in my code.

location = 'Times Square, New York, NY'
when = '2023-01-01 00:00'

First, let’s get the latitude and longitude of our location. We will use the geopy library to do that.

locator = Nominatim(user_agent='myGeocoder')
location = locator.geocode(location)
lat, long = location.latitude, location.longitude

Second, using tzwhere and pytz libraries, get the timezone of our location based on the latitude and longitude and date-time.

# convert date string into datetime object
dt = datetime.strptime(when, '%Y-%m-%d %H:%M')
# define datetime and convert to utc based on our timezone
timezone_str = tzwhere.tzwhere().tzNameAt(lat, long)
local = timezone(timezone_str)
# get UTC from local timezone and datetime
local_dt = local.localize(dt, is_dst=None)
utc_dt = local_dt.astimezone(utc)

Third, using the skyfield library, define the location and time of our observation and define the position of our observer. We use world geodetic data in order to generate our observer object.

# find location of earth and sun and set the observer position
sun = eph[‘sun’]
earth = eph[‘earth’]
# define observation time from our UTC datetime
ts = load.timescale()
t = ts.from_datetime(utc_dt)
# define an observer using the world geodetic system data
observer = wgs84.latlon(latitude_degrees=lat, longitude_degrees=long).at(t)
position = observer.from_altaz(alt_degrees=90, az_degrees=0)

*Note: I think we don’t actually need to define the position here because observer defaults to this view but I include it here in case you would like to change the angle and direction of your observation.

Calculate Star Positions

At this point, we finally have our location. We can now use it to calculate the location of stars as they would appear in the sky.

To do so, we need to give skyfield library something to observe. We will create a fake star right above our heads (at latitude of 90 degrees) and center our observation on it. The radec method generates the right ascension and declination from the International Celestial Reference System (ICRS) and we use it to define where the star is.

ra, dec, distance = observer.radec()
center_object = Star(ra=ra, dec=dec)

Now we observe the position of our star in space from Earth’s location in space and build a stereographic projection of the sky. A stereographic projection essentially converts a sphere into a plane field.

center = earth.at(t).observe(center_object)
projection = build_stereographic_projection(center)

Finally, let’s observe all of stars from Earth and project them onto our sky map as well. This will calculate their positions in terms of x and y on our 2d plane.

star_positions = earth.at(t).observe(Star.from_dataframe(stars))
stars[‘x’], stars[‘y’] = projection(star_positions)

Build the Star Chart

I will use matplotlib to build the chart but you can use any number of alternatives (like pillow, for example).

Calculating star size based on its apparent magnitude

To make the star chart more representative of the actual sky, we need to make the stars variable in brightness and to do so in our rudimentary 2d plane, we can just vary their size. Stars have varying apparent magnitudes which define how bright they are compared the brightest star in our sky. The greater the magnitude, the less bright the star appears to us.

A formula for apparent magnitude, based on observed relative brightness of stars is

Where Fx represents observed irradiance using spectral filter and Fx0 is its reference flux. We can change this formula to suit our need.

Let m_i represent the apparent magnitude of our stars (which we already have in the stars dataset). Let S_i represent the size of star i in our chart and S_0 represent the maximum size of a star in our chart. We arbitrarily select S_0 depending how big we want to go. Then,

And solving for S_i, we get

Now we can use this formula to define the size of our stars.

Plotting the Star Chart

I define the size of my star chart as 10x10 (inches) and set the maximum star size as 100.

We probably don’t want to see ALL of the stars so we should limit which ones we actually produce using limiting_magnitude parameter. Let’s assume it’s 10 (which includes a lot of stars and definitely not realistic for NYC). We then filter our stars data to only include stars which are as bright or brighter than that with magnitude 10.

chart_size = 10
max_star_size = 100
limiting_magnitude = 10
bright_stars = (stars.magnitude <= limiting_magnitude)
magnitude = stars[‘magnitude’][bright_stars]

I want to generate a dark sky with white stars so I will create a matplotlib figure and add a dark circle with radius of 1, centered at (0, 0) as a background for my stars.

fig, ax = plt.subplots(figsize=(chart_size, chart_size))
 
border = plt.Circle((0, 0), 1, color=’navy’, fill=True)
ax.add_patch(border)

Using my formula above, I calculate my star marker size:

marker_size = max_star_size * 10 ** (magnitude / -2.5)

I create a scatterplot with my stars, using their x and y location, marker size which represents brightness. I use circle marker type without any outline and I set zorder=2 to ensure that my stars appear on top of the navy circle.

ax.scatter(stars[‘x’][bright_stars], stars[‘y’][bright_stars],
 s=marker_size, color=’white’, marker=’.’, linewidths=0, 
 zorder=2)

I will also clip the horizon (this will remove everything outside of my circle). It’s not so helpful for me since my stars are white but would be if I put black stars on white background instead.

horizon = Circle((0, 0), radius=1, transform=ax.transData)
for col in ax.collections:
    col.set_clip_path(horizon)

Finally, I remove axes, set axis limits, and plot the chart:

ax.set_xlim(-1, 1)
ax.set_ylim(-1, 1)
plt.axis(‘off’)
plt.show()

The Final Product

Voila! Now you should have a star map like this one:

The final star map

You can clearly see some cool stuff in the star map above: not the Orion’s belt around 5 o’clock and the big dipper at 10:30!

You can check out the source code in a jupyter notebook here.

How can we make this even better?

• We can add altitude and azimuth lines
• We can add a fisheye effect in order to make the sky more realistic from the observer point of view
• We can add constellations
• We can add labels for stars and constellations
• We can add the moon and the sun
• We can add comets
• We can observe the sky at different angles
• We can place it on the light-to-dark gradient aligned to east and west positions relative to the direction of observation and time
• We can cut it out of wood using a laser cutting machine or cut it out in vinyl using a vinyl cutting machine
• We can put it on a t-shirt or make an art print
• We can make a gif which shows star movement over a 24 hour period

These ideas are trivial and are left as an exercise to the reader.