import astropy.coordinates as coord
import astropy.table as at
import astropy.units as u
import gala.coordinates as gc
import healpy as hp
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
from healpy.projaxes import HpxMollweideAxes
import pyia

%matplotlib inline

Load DESI MWS DR1 data and positional cross-match to Gaia DR3:

desi = at.Table.read("~/data/DESI/DR1/mwsall-pix-iron.fits")
gaia = at.Table.read("~/data/DESI/DR1/desi_dr1_gaia_dr3_xm.fits")
len(desi), len(gaia)

WARNING: hdu= was not specified but multiple tables are present, reading in first available table (hdu=1) [astropy.io.fits.connect]
WARNING: UnitsWarning: 'log(cm.s**-2)' did not parse as fits unit: 'log' is not a recognized function If this is meant to be a custom unit, define it with 'u.def_unit'. To have it recognized inside a file reader or other code, enable it with 'u.add_enabled_units'. For details, see https://docs.astropy.org/en/latest/units/combining_and_defining.html [astropy.units.core]
WARNING: UnitsWarning: 'log(cm.s**-2)' did not parse as fits unit: 'log' is not a recognized function If this is meant to be a custom unit, define it with 'u.def_unit'. To have it recognized inside a file reader or other code, enable it with 'u.add_enabled_units'. For details, see https://docs.astropy.org/en/latest/units/combining_and_defining.html [astropy.units.core]
WARNING: UnitsWarning: 'log(cm.s**-2)' did not parse as fits unit: 'log' is not a recognized function If this is meant to be a custom unit, define it with 'u.def_unit'. To have it recognized inside a file reader or other code, enable it with 'u.add_enabled_units'. For details, see https://docs.astropy.org/en/latest/units/combining_and_defining.html [astropy.units.core]
WARNING: UnitsWarning: 'log(cm.s**-2)' did not parse as fits unit: 'log' is not a recognized function If this is meant to be a custom unit, define it with 'u.def_unit'. To have it recognized inside a file reader or other code, enable it with 'u.add_enabled_units'. For details, see https://docs.astropy.org/en/latest/units/combining_and_defining.html [astropy.units.core]

(6372607, 5381083)

There are some duplicates in the cross-match -- either because a source was observed twice in DESI, or because there are multiple Gaia sources within the DESI fiber footprint:

unq, counts = np.unique(gaia["targetid"], return_counts=True)
duplicate_targetids = unq[counts > 1]
(counts > 1).sum()

np.int64(212471)

For now, we'll ignore the duplicates:

joined = at.join(
    desi[~np.isin(desi["TARGETID"], duplicate_targetids)],
    gaia,
    keys_left="TARGETID",
    keys_right="targetid",
    join_type="inner",
)
len(joined), (counts == 1).sum()

(4940396, np.int64(4940376))

Since we are working with Gaia data, we can use the pyia package to help:

g = pyia.GaiaData(
    joined, radial_velocity_colname="VRAD", radial_velocity_error_colname="VRAD_ERR"
)

Sky plots:

First make some visualizations of the DESI data on the sky (using healpy):

c = g.get_skycoord(distance=False)

def apw_healpix_plot(
    c: coord.SkyCoord,
    nside: int = 128,
    figsize: tuple[int, int] = (12, 6),
    projmap_kwargs: dict | None = None,
    graticule: bool = True,
):
    hpix = hp.ang2pix(nside, c.data.lon.degree, c.data.lat.degree, lonlat=True)
    hpix_counts = np.bincount(hpix, minlength=hp.nside2npix(nside))

    projmap_kwargs = projmap_kwargs or dict()
    projmap_kwargs.setdefault("vmin", 1.0)
    projmap_kwargs.setdefault("vmax", np.nanmax(hpix_counts))
    projmap_kwargs.setdefault("cmap", "magma_r")
    projmap_kwargs.setdefault("norm", mpl.colors.LogNorm())
    projmap_kwargs.setdefault("badcolor", "w")

    fig = plt.figure(figsize=figsize)
    proj_ax = HpxMollweideAxes(fig, 111)
    proj_ax.projmap(hpix_counts, **projmap_kwargs)

    if graticule:
        proj_ax.graticule(30.0, 30.0, color="tab:blue", alpha=0.75)
    fig.add_axes(proj_ax)

    return fig

fig = apw_healpix_plot(
    c,
    projmap_kwargs=dict(
        rot=[180, 0, 0],
    ),
)
fig.suptitle("ICRS centered on RA=180º");

Now we can look at the sky coverage in different coordinate systems, like Galactic, and a few different stellar stream coordinate frames:

for frame in [
    coord.Galactic(),
    gc.SagittariusVasiliev21(),
    gc.GD1Koposov10(),
    gc.Pal5PriceWhelan18()
]:
    c_fr = c.transform_to(frame)
    fig = apw_healpix_plot(c_fr)
    fig.suptitle(frame.__class__.__name__)

Most of the DESI data does not cover my favorite streams, but it seems like Pal 5 may be in the footprint! LEt's look for some serendipitously-observed stars?

Pal 5

To get a sense of the footprint, I load data for RR Lyrae stars in the Pal 5 stream (from Price-Whelan et al. 2019):

rrl_pal5_g = pyia.GaiaData("~/projects/pal5-rrl/data/RRL-with-prob.fits")
rrl_pal5_g = rrl_pal5_g[
    (rrl_pal5_g.member_prob > 0.8) & (rrl_pal5_g.inside_stream_track)
]
c_rrl_pal5 = rrl_pal5_g.get_skycoord(distance=False).transform_to(
    gc.Pal5PriceWhelan18()
)
len(rrl_pal5_g)

24

c_pal5 = c.transform_to(gc.Pal5PriceWhelan18())

Some rough selections to get stream stars:

mask = (
    (np.abs(c_pal5.phi1.wrap_at(180 * u.deg)) < 25 * u.deg)
    & (c_pal5.phi2 > -2 * u.deg)
    & (c_pal5.phi2 < 5 * u.deg)
    & (np.abs(c.pm_ra_cosdec.value - -2.75) < 1.5)
    & (np.abs(c.pm_dec.value - -2.69) < 1.5)
    & (g.FEH < -0.5)
)

plt.hist2d(
    g.bp_rp[mask].value,
    g.phot_g_mean_mag[mask].value,
    bins=(np.linspace(-0.5, 3.5, 128), np.linspace(12, 21, 128))
);

c_pal5_masked = c_pal5[mask]
len(c_pal5_masked)

8705

plt.figure(figsize=(10, 5))
plt.plot(c_pal5_masked.phi1.degree, c_pal5_masked.phi2.degree, ls="none", alpha=0.5)
plt.scatter(c_rrl_pal5.phi1.degree, c_rrl_pal5.phi2.degree, color="tab:blue")
plt.xlim(-25, 25)
plt.ylim(-2, 5)

(-2.0, 5.0)

plt.figure(figsize=(10, 5))
plt.plot(
    c_pal5_masked.phi1.degree, c_pal5_masked.radial_velocity.value, ls="none", alpha=0.5
)
plt.ylim(-150, 50)

(-150.0, 50.0)