# Licensed under a 3-clause BSD style license - see LICENSE.rst
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import corner
import arviz_plots
import astropy.units as u
from astropy.visualization import hist
from pytelpoint.stats import psd, skyrms
from pytelpoint.fitting import best_fit_pars
__all__ = [
"pointing_histogram",
"pointing_residuals",
"pointing_sky",
"pointing_azel_resid",
"plot_corner",
"plot_posterior",
]
[docs]
def pointing_histogram(coo_ref, coo_meas, bins="freedman", style="ggplot"):
"""
Plot histogram of separations between reference coordinates, coo_ref, and measured coordinates, coo_meas.
Parameters
----------
coo_ref : `~astropy.coordinates.SkyCoord` instance
Reference coordinates
coo_meas : `~astropy.coordinates.SkyCoord` instance
Measured coordinates
bins : str (default: freedman)
Name of algorithm to use to calculate optimal bin sizes. See `~astropy.visualization.hist` for options
and descriptions.
style : str (default: ggplot)
Matplotlib style to apply to the plot
Returns
-------
fig : `matplotlib.figure.Figure` instance
Figure object containing the histogram plot.
"""
seps = coo_ref.separation(coo_meas)
with plt.style.context(style, {"xtick.labelsize": 18, "ytick.labelsize": 18}):
matplotlib.rcParams.update({"font.size": 15})
fig, ax = plt.subplots(figsize=[9, 6])
hist(
seps.to(u.arcsec).value, bins=bins, ax=ax, histtype="stepfilled", alpha=0.6
)
ax.set_ylabel("N")
ax.set_xlabel("Pointing Error (arcsec)")
med = np.median(seps.to(u.arcsec))
rms = skyrms(coo_ref, coo_meas)
skypsd = psd(coo_ref, coo_meas)
ax.vlines(
x=med.value,
ymin=0,
ymax=1,
transform=ax.get_xaxis_transform(),
color="gray",
alpha=0.5,
ls="-",
label=f'Median: {med.value:.2f}"',
)
ax.vlines(
x=rms.value,
ymin=0,
ymax=1,
transform=ax.get_xaxis_transform(),
color="gray",
alpha=0.5,
ls="--",
label=f'RMS: {rms.value:.2f}"',
)
ax.vlines(
x=skypsd.value,
ymin=0,
ymax=1,
transform=ax.get_xaxis_transform(),
color="gray",
alpha=0.5,
ls=":",
label=f'PSD: {skypsd.value:.2f}"',
)
ax.legend()
plt.tight_layout()
return fig
[docs]
def pointing_residuals(coo_ref, coo_meas, circle_size=1.0, style="ggplot"):
"""
Plot of pointing residuals in az/el space
Parameters
----------
coo_ref : `~astropy.coordinates.SkyCoord` instance
Reference coordinates
coo_mes : `~astropy.coordinates.SkyCoord` instance
Measured coordinates
circle_size : float (default: 1.)
Size of reference circle to plot with residuals. Set to None
for no circle.
style : str (default: ggplot)
Matplotlib style to apply to the plot
Returns
-------
fig : `matplotlib.figure.Figure` instance
Figure object containing the residuals plot.
"""
az_res, el_res = coo_meas.spherical_offsets_to(coo_ref)
with plt.style.context(style, {"xtick.labelsize": 18, "ytick.labelsize": 18}):
matplotlib.rcParams.update({"font.size": 16})
fig, ax = plt.subplots(figsize=[6, 6])
ax.set_aspect("equal")
ax.scatter(az_res.to(u.arcsec), el_res.to(u.arcsec))
if circle_size is not None:
ax_size = 5 * circle_size
ax.set_xlim([-ax_size, ax_size])
ax.set_ylim([-ax_size, ax_size])
ax.set_xlabel(r"$\Delta$A (arcsec)")
ax.set_ylabel(r"$\Delta$E (arcsec)")
ax.set_title("Azimuth-Elevation Residuals")
if circle_size is not None:
c1 = matplotlib.patches.Circle(
(0, 0),
circle_size,
ec="black",
lw=4,
fill=False,
alpha=0.4,
label=f'{circle_size}"',
)
ax.add_patch(c1)
ax.legend()
plt.tight_layout()
return fig
[docs]
def pointing_sky(coo_ref, coo_meas, style="ggplot"):
"""
Plot of pointing errors as a function of sky position
Parameters
----------
coo_ref : `~astropy.coordinates.SkyCoord` instance
Reference coordinates
coo_mes : `~astropy.coordinates.SkyCoord` instance
Measured coordinates
style : str (default: ggplot)
Matplotlib style to apply to the plot
Returns
-------
fig : `matplotlib.figure.Figure` instance
Figure object containing the pointing errors plot.
"""
az_res, el_res = coo_meas.spherical_offsets_to(coo_ref)
with plt.style.context(style, {"xtick.labelsize": 16, "ytick.labelsize": 16}):
matplotlib.rcParams.update({"font.size": 14})
x = coo_ref.az
y = 90 * u.degree - coo_ref.alt # use zenith angle here as a trick
uu = (az_res).to(u.arcsec).value
vv = (-1 * el_res).to(u.arcsec).value
fig, ax = plt.subplots(subplot_kw={"projection": "polar"}, figsize=[8, 8])
qq = ax.quiver(
x.to(u.radian).value,
y.value,
uu,
vv,
np.sqrt(uu**2 + vv**2),
scale_units="y",
angles="xy",
pivot="tip",
color="red",
)
ax.set_rmax(90)
ticks = [0, 15, 30, 45, 60, 75, 90]
ax.set_rticks(ticks)
ax.set_rlim(0, 90)
ax.set_yticklabels(
[
r"$90^{\circ}$",
r"$75^{\circ}$",
r"$60^{\circ}$",
r"$45^{\circ}$",
r"$30^{\circ}$",
r"$15^{\circ}$",
r"El = $0^{\circ}$",
]
)
ax.set_theta_zero_location("N")
ax.set_theta_direction(-1)
ax.scatter(x.to(u.radian).value, y.value)
cbar = plt.colorbar(qq, shrink=0.7)
cbar.set_label("arcsec")
plt.tight_layout()
return fig
[docs]
def pointing_azel_resid(coo_ref, coo_meas, style="ggplot"):
"""
Plot of az/el pointing residuals as a function of az/el
Parameters
----------
coo_ref : `~astropy.coordinates.SkyCoord` instance
Reference coordinates
coo_mes : `~astropy.coordinates.SkyCoord` instance
Measured coordinates
style : str (default: ggplot)
Matplotlib style to apply to the plot
Returns
-------
fig : `matplotlib.figure.Figure` instance
Figure object containing the az/el residuals plot.
"""
az_res, el_res = coo_meas.spherical_offsets_to(coo_ref)
az_res = az_res.to(u.arcsec)
el_res = el_res.to(u.arcsec)
az = coo_ref.az
el = coo_ref.alt
azel_max = np.max(
[
5,
np.abs(az_res).max().to(u.arcsec).value,
np.abs(el_res).max().to(u.arcsec).value,
]
)
with plt.style.context(style, {"xtick.labelsize": 18, "ytick.labelsize": 18}):
matplotlib.rcParams.update({"font.size": 16})
fig, axs = plt.subplots(2, 2, figsize=(12, 8), sharex="col", sharey="row")
axs[0, 0].set_ylim(-azel_max, azel_max)
axs[1, 0].set_xlim(0, 360)
axs[1, 0].set_ylim(-azel_max, azel_max)
axs[1, 1].set_xlim(0, 90)
axs[0, 0].scatter(az, az_res)
axs[1, 0].scatter(az, el_res)
axs[0, 1].scatter(el, az_res)
axs[1, 1].scatter(el, el_res)
axs[0, 0].set_ylabel(r"$\Delta A$ (arcsec)")
axs[1, 0].set_ylabel(r"$\Delta E$ (arcsec)")
axs[1, 0].set_xlabel("Azimuth")
axs[1, 0].set_xticks([0, 90, 180, 270, 360])
axs[1, 0].set_xticklabels(["N", "E", "S", "W", "N"])
axs[1, 1].set_xlabel("Elevation")
axs[1, 1].set_xticks([0, 15, 30, 45, 60, 75, 90])
axs[1, 1].set_xticklabels(
[f"{i}" + r"$^{\circ}$" for i in [0, 15, 30, 45, 60, 75, 90]]
)
plt.tight_layout()
return fig
[docs]
def plot_corner(
idata,
quantiles=None,
truths=None,
title_kwargs={"fontsize": 18},
label_kwargs={"fontsize": 16},
):
"""
Make corner plot from outputs of a pymc az/el fit
Parameters
----------
idata : object
Any object that can be converted to an `~arviz.InferenceData` object. Refer to documentation of
arviz.convert_to_dataset for details.
quantiles : list of float (default: None -> [0.16, 0.5, 0.84])
Quantiles to overlay on each histogram plot
truths : dict
Dict of reference parameters to overlay on plots. Must contain the following keys:
ia, ie, an, aw, ca, npae, tf, tx, el_sigma, and az_sigma.
Values set to None won't be displayed. Default is to not display any.
Returns
-------
fig : `matplotlib.figure.Figure` instance
Figure object containing the corner plot.
"""
if truths is None:
truths = {}
if quantiles is None:
quantiles = [0.16, 0.5, 0.84]
pars = best_fit_pars(idata)
labels = []
for p in pars.keys():
if p not in truths:
truths[p] = None
if "sigma" in p:
ax, _ = p.split("_")
labels.append("$\\sigma_{" + ax.upper() + "}$")
else:
labels.append(p.upper())
fig = corner.corner(
idata,
labels=labels,
quantiles=quantiles,
truths=truths,
show_titles=True,
title_kwargs=title_kwargs,
label_kwargs=label_kwargs,
)
return fig
[docs]
def plot_posterior(idata):
"""
Make posterior probability distributions plot from a pymc fit
Parameters
----------
idata : object
Any object that can be converted to an `~arviz.InferenceData` object.
Refer to documentation of arviz.convert_to_dataset for details.
Returns
-------
pc : `~arviz_plots.PlotCollection`
PlotCollection containing the posterior distribution plots.
"""
pc = arviz_plots.plot_dist(idata)
return pc