Possible bug in the plot_sky function

[gnthibault]

Dear all,

First I would like to congratulate the developers of this awesome project. I was in the process of developing a similar thing, but I was very happy to see that someone already have implemented many useful features here. Long life to this project !

I took a look at the plot_sky function, but there is a line of code I am not sure I understand:

line 134:

altitude = (91 * u.deg - observer.altaz(time, target).alt) * (1/u.deg) Why are you using 91 degrees instead of 90 ? An object at 0.5 elevation would be remapped to 90.5, hence below the horizon on the drawing ?

Please correct me if I am wrong

Thank you in advance for your help

[bmorris3]

I'm so sorry we let this issue slip through the cracks for so many months. Thank you for your kind words 🎉

I can't remember how this number got to 91, @eteq – do you remember what edge case we were avoiding?

[bmorris3]

Hi @gnthibault,

I dug into this a bit, and it's a visualization convention that we adopted. We did this because for plots with project='polar', the plt.plot(phi, theta) command would, of course, plot theta as a function of phi. That's what most people expect of a polar plotter, so that's what matplotlib does.

Now let's suppose we call some function plot_sky_map(az, alt). In a plot of the sky, we usually visualize alt=0, on the horizon, which we represent with a point on the edge of the circular plot. alt=90 deg would be at the zenith, which we normally represent with the center of the circular plot. So a first guess for the transformation necessary to plot the sky "astronomer" style is to do plt.plot(az, (90*u.deg - alt).value). This mostly works, unless targets fall right near the horizon. If you want to plot targets on the horizon, you have to be careful, and I think that's where 91 comes in.

I coded up a toy example that demonstrates that (a) the precise value 91 doesn't matter much, and (b) under certain conditions it provides better behavior than 90. See below!

import matplotlib.pyplot as plt
import numpy as np
import astropy.units as u

from astroplan import Observer
from astropy.coordinates import SkyCoord
from astropy.time import Time

# cycle through the first several stars in Orion:
greek_letters = 'alf bet del gam eps zet eta'.split()
constellation = 'Ori'

# This time and observer puts the targets near the horizon.
# Use an hour or two earlier to put the stars above horizon.
observer = Observer.at_site("APO")
time = Time("2022-10-14 17:40")

targets = []

for name in [f"{letter} {constellation}" for letter in greek_letters]:
    targets.append(SkyCoord.from_name(name))
    
targets = SkyCoord(targets)

# Compute the azimuth and altitude in units of radians
azimuth = observer.altaz(time, targets).az.rad
altitude = observer.altaz(time, targets).alt.rad

# Try different constants, show that the result works either way:
for try_constant in [90, 91]:

    rescaled_altitude = (np.radians(try_constant) - altitude)

    fig = plt.figure(figsize=(10, 4))
    fig.suptitle(f'constant = ${try_constant}^\\circ$')
    ax_linear = fig.add_subplot(121)
    ax_polar = fig.add_subplot(122, projection='polar')

    for axis in [ax_linear, ax_polar]:
        axis.scatter(azimuth, rescaled_altitude)

    # Set radial limits:
    ax_polar.set_rlim(np.radians([try_constant - 90, try_constant]))
    ax_linear.axhline(np.radians(try_constant))
    ax_linear.set(
        xlabel="az [rad]", ylabel="rescaled alt [rad]"
    )

Hope this helps. Thanks!