Python – How to Plot Star Maps Using the Equatorial Coordinate System with Matplotlib, Astropy, and Skyfield

astropymatplotlibpythonskyfield

I'm trying to generate star maps with the equatorial coordinates system (RAJ2000 and DEJ2000). However, I only get a grid system where meridians and parallels are in parallel, while parallels should be curved and meridians should converge to the north celestial pole and the south ceestial pole.

I'm using some Python modules: matplotlib, skyfield (for the stereographic projection), astroquery (so I can target any object in the deep space) and astropy.

This is my code:

#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""Generate a skymap with equatorial grid"""

import numpy as np
from matplotlib import pyplot as plt
from matplotlib.collections import LineCollection
from skyfield.api import Star, load
from skyfield.data import hipparcos, stellarium
from skyfield.projections import build_stereographic_projection
from astroquery.simbad import Simbad
from astropy.coordinates import SkyCoord
import astropy.units as u
from astropy.wcs import WCS
from astropy.visualization.wcsaxes import WCSAxes

# Design
plt.style.use("dark_background")
plt.rcParams['font.family'] = 'serif'
plt.rcParams['font.serif'] = ['Times New Roman']

# Query object from Simbad
OBJECT = "Alioth"
FOV = 30.0
MAG = 6.5

TABLE = Simbad.query_object(OBJECT)
RA = TABLE['RA'][0]
DEC = TABLE['DEC'][0]
COORD = SkyCoord(f"{RA} {DEC}", unit=(u.hourangle, u.deg), frame='fk5')

print("RA is", RA)
print("DEC is", DEC)

ts = load.timescale()
t = ts.now()

# An ephemeris from the JPL provides Sun and Earth positions.
eph = load('de421.bsp')
earth = eph['earth']

# Load constellation outlines from Stellarium
url = ('https://raw.githubusercontent.com/Stellarium/stellarium/master'
       '/skycultures/modern_st/constellationship.fab')

with load.open(url) as f:
    constellations = stellarium.parse_constellations(f)

edges = [edge for name, edges in constellations for edge in edges]
edges_star1 = [star1 for star1, star2 in edges]
edges_star2 = [star2 for star1, star2 in edges]

# The Hipparcos mission provides our star catalog.
with load.open(hipparcos.URL) as f:
    stars = hipparcos.load_dataframe(f)

# Center the chart on the specified object's position.
center = earth.at(t).observe(Star(ra_hours=COORD.ra.hour, dec_degrees=COORD.dec.degree))
projection = build_stereographic_projection(center)

# Compute the x and y coordinates that each star will have on the plot.
star_positions = earth.at(t).observe(Star.from_dataframe(stars))
stars['x'], stars['y'] = projection(star_positions)

# Create a True/False mask marking the stars bright enough to be included in our plot.
bright_stars = (stars.magnitude <= MAG)
magnitude = stars['magnitude'][bright_stars]
marker_size = (0.5 + MAG - magnitude) ** 2.0

# The constellation lines will each begin at the x,y of one star and end at the x,y of another.
xy1 = stars[['x', 'y']].loc[edges_star1].values
xy2 = stars[['x', 'y']].loc[edges_star2].values
lines_xy = np.rollaxis(np.array([xy1, xy2]), 1)

# Define the limit for the plotting area
angle = np.deg2rad(FOV / 2.0)
limit = np.tan(angle)  # Calculate limit based on the field of view

# Build the figure with WCS axes
fig = plt.figure(figsize=[6, 6])
wcs = WCS(naxis=2)
wcs.wcs.crpix = [1, 1]
wcs.wcs.cdelt = np.array([-FOV / 360, FOV / 360])
wcs.wcs.crval = [COORD.ra.deg, COORD.dec.deg]
wcs.wcs.ctype = ["RA---STG", "DEC--STG"]

ax = fig.add_subplot(111, projection=wcs)

# Draw the constellation lines
ax.add_collection(LineCollection(lines_xy, colors='#ff7f2a', linewidths=1, linestyle='-'))

# Draw the stars
ax.scatter(stars['x'][bright_stars], stars['y'][bright_stars],
           s=marker_size, color='white', zorder=2)

ax.scatter(RA, DEC, marker='*', color='red', zorder=3)

angle = np.pi - FOV / 360.0 * np.pi
limit = np.sin(angle) / (1.0 - np.cos(angle))

# Set plot limits
ax.set_xlim(-limit, limit)
ax.set_ylim(-limit, limit)
ax.set_aspect('equal')

# Add RA/Dec grid lines
ax.coords.grid(True, color='white', linestyle='dotted')

# Set the coordinate grid
ax.coords[0].set_axislabel('RA (hours)')
ax.coords[1].set_axislabel('Dec (degrees)')
ax.coords[0].set_major_formatter('hh:mm:ss')
ax.coords[1].set_major_formatter('dd:mm:ss')

# Title
ax.set_title(f'Sky map centered on {OBJECT}', color='white', y=1.04)

# Save the image
FILE = "chart.png"
plt.savefig(FILE, dpi=100, facecolor='#1a1a1a')

And this is the resulting image:

As you can see, the grid (parallels and meridians) are totally parallel. However, my goal is to achieve this grid:

In this case, I got the right WCS from a FITS image in the DSS survey. It's automatic. However, for plotting star maps I need to create a simulation of that, working fine with the labels and the coordinates system, not as a background image or something else.

Best Answer

I can't replicate your code, to show where exactly you need corrections. Matplotlib allows you to use polar coordinates,for this you need ti apply axes and grids like in those docs: matplotlib docs

If you can simplify code without using files (with just list of numbers or strings, idk), I can try to fix it

The following is a simple example of how you can do this:

    #!/usr/bin/env python
# -*- coding: utf-8 -*-
"""Generate a skymap with equatorial grid"""

import numpy as np
from matplotlib import pyplot as plt
from matplotlib.collections import LineCollection
from skyfield.api import Star, load
from skyfield.data import hipparcos, stellarium
from skyfield.projections import build_stereographic_projection
from astroquery.simbad import Simbad
from astropy.coordinates import SkyCoord
import astropy.units as u
from astropy.wcs import WCS
from astropy.visualization.wcsaxes import WCSAxes
import numpy as np

from matplotlib.projections import PolarAxes
from matplotlib.transforms import Affine2D
from mpl_toolkits.axisartist import Axes, HostAxes, angle_helper
from mpl_toolkits.axisartist.grid_helper_curvelinear import \
    GridHelperCurveLinear

# Design
plt.style.use("dark_background")
plt.rcParams['font.family'] = 'serif'
plt.rcParams['font.serif'] = ['Times New Roman']

# Query object from Simbad
OBJECT = "Polaris"
FOV = 30.0
MAG = 6.5

TABLE = Simbad.query_object(OBJECT)
RA = TABLE['RA'][0]
DEC = TABLE['DEC'][0]
COORD = SkyCoord(f"{RA} {DEC}", unit=(u.hourangle, u.deg), frame='fk5')

print("RA is", RA)
print("DEC is", DEC)

ts = load.timescale()
t = ts.now()

# An ephemeris from the JPL provides Sun and Earth positions.
eph = load('de421.bsp')
earth = eph['earth']

# Load constellation outlines from Stellarium
url = ('https://raw.githubusercontent.com/Stellarium/stellarium/master'
       '/skycultures/modern_st/constellationship.fab')

with load.open(url) as f:
    constellations = stellarium.parse_constellations(f)

edges = [edge for name, edges in constellations for edge in edges]
edges_star1 = [star1 for star1, star2 in edges]
edges_star2 = [star2 for star1, star2 in edges]

# The Hipparcos mission provides our star catalog.
with load.open(hipparcos.URL) as f:
    stars = hipparcos.load_dataframe(f)

# Center the chart on the specified object's position.
center = earth.at(t).observe(Star(ra_hours=COORD.ra.hour, dec_degrees=COORD.dec.degree))
projection = build_stereographic_projection(center)

# Compute the x and y coordinates that each star will have on the plot.
star_positions = earth.at(t).observe(Star.from_dataframe(stars))
stars['x'], stars['y'] = projection(star_positions)

# Create a True/False mask marking the stars bright enough to be included in our plot.
bright_stars = (stars.magnitude <= MAG)
magnitude = stars['magnitude'][bright_stars]
marker_size = (0.5 + MAG - magnitude) ** 2.0

# The constellation lines will each begin at the x,y of one star and end at the x,y of another.
xy1 = stars[['x', 'y']].loc[edges_star1].values
xy2 = stars[['x', 'y']].loc[edges_star2].values
lines_xy = np.rollaxis(np.array([xy1, xy2]), 1)

# Define the limit for the plotting area
angle = np.deg2rad(FOV / 2.0)
limit = np.tan(angle)  # Calculate limit based on the field of view

# Build the figure with WCS axes
fig = plt.figure(figsize=[6, 6])
wcs = WCS(naxis=2)
wcs.wcs.crpix = [1, 1]
wcs.wcs.cdelt = np.array([-FOV / 360, FOV / 360])
wcs.wcs.crval = [COORD.ra.deg, COORD.dec.deg]
wcs.wcs.ctype = ["RA---STG", "DEC--STG"]

# <- Start point code implementetion from docs
tr = Affine2D().scale(np.pi/180, 1) + PolarAxes.PolarTransform()
grid_locator1 = angle_helper.LocatorDMS(12)

grid_helper = GridHelperCurveLinear(
    tr, grid_locator1=grid_locator1)

ax = fig.add_subplot(
    1, 1, 1, axes_class=HostAxes, grid_helper=grid_helper)
ax.axis["right"].major_ticklabels.set_visible(True)
ax.axis["top"].major_ticklabels.set_visible(True)
ax.axis["right"].get_helper().nth_coord_ticks = 0
ax.axis["bottom"].get_helper().nth_coord_ticks = 1

ax.set_aspect(1)
ax.grid(True, zorder=0)
ax.grid(True, color='white', linestyle='dotted')
# <- End point code implementetion from docs

# Draw the constellation lines
ax.add_collection(LineCollection(lines_xy, colors='#ff7f2a', linewidths=1, linestyle='-'))

# Draw the stars
ax.scatter(stars['x'][bright_stars], stars['y'][bright_stars],
           s=marker_size, color='white', zorder=2)

ax.scatter(RA, DEC, marker='*', color='red', zorder=3)

angle = np.pi - FOV / 360.0 * np.pi
limit = np.sin(angle) / (1.0 - np.cos(angle))

# Set plot limits
# Also some changes below
ax.set_xlim(-limit, limit)
ax.set_ylim(-limit, limit)
ax.set_aspect('equal')

# Add RA/Dec grid lines
ax.grid(True, color='white', linestyle='dotted')

# Set the coordinate grid
ax.set_xlabel('RA (hours)')
ax.set_ylabel('Dec (degrees)')
ax.xaxis.set_major_formatter('hh:mm:ss')
ax.yaxis.set_major_formatter('dd:mm:ss')

# Title
ax.set_title(f'Sky map centered on {OBJECT}', color='white', y=1.04)

# Save the image
FILE = "chart.png"
plt.savefig(FILE, dpi=100, facecolor='#1a1a1a')

Maybe I missed something, but it seems ok

Note:

A Python dict, semantically used for keyword argument passing, is arbitrarily ordered. However, in Python 3.6+, keyword arguments are guaranteed to remember insertion order. "The order of elements in **kwargs now corresponds to the order in which keyword arguments were passed to the function." - What’s New In Python 3.6. In fact, all dicts in CPython 3.6 will remember insertion order as an implementation detail, and this becomes standard in Python 3.7.

Best Answer

Related Solutions

Python Subprocess – How to Execute a Program or Call a System Command

Python ** and * Parameters – What Do They Do?

Note:

Related Question