Y4 Y5 MWS program selection

[segasai]

Hi,

This is the starting point to iterate over the integration of the selection of GD1 targets into desitarget following the discussion yesterday. (@geordie666 , @sazabi4 )

The code below is based on my notebook and is self-contained, in the sense that I integrated several dependencies there (i.e. various rotation, proper motion transformations from github.com/segasai/astrolibpy, specifically ( i.e.https://github.com/segasai/astrolibpy/blob/master/my_utils/rotate_pm.py https://github.com/segasai/astrolibpy/blob/master/my_utils/sphere_rotate.py ).

The main selection code is in the main section.
The reading of the data is in read_data() and that's what need to be replaced by the DECals X Gaia DR3 match from desitarget.

The only non-standard dependency that I cannot easily get away from is the gaia_zeropoint
https://gitlab.com/icc-ub/public/gaiadr3_zeropoint
which is necessary to fix the offset in Gaia parallaxes.

I can probably create a PR based on this code, provided there are no major objections to how things are structured now.

import healpy
import datetime
import os
import scipy.interpolate
import numpy as np
import astropy.table as atpy

import astropy.coordinates as acoo
import astropy.units as auni
from zero_point import zero_point as gaia_zpt

# Use v4.0 defaults
GCPARAMS = acoo.galactocentric_frame_defaults.get_from_registry(
    "v4.0")['parameters']

kms = auni.km / auni.s
masyr = auni.mas / auni.year

gaia_zpt.zpt.load_tables()

def cosd(x):
    return np.cos(np.deg2rad(x))


def sind(x):
    return np.sin(np.deg2rad(x))

def betw(x, x1, x2):
    return (x >= x1) & (x < x2)


def torect(ra, dec):
    x = np.cos(ra / 57.295779513082323) * np.cos(dec / 57.295779513082323)
    y = np.sin(ra / 57.295779513082323) * np.cos(dec / 57.295779513082323)
    z = np.sin(dec / 57.295779513082323)
    return x, y, z


def fromrect(x, y, z):
    ra = np.arctan2(y, x) * 57.295779513082323
    dec = 57.295779513082323 * np.arctan2(z, np.sqrt(x**2 + y**2))
    return ra, dec


def rotation_matrix(rapol, decpol, ra0):
    """
    Return the rotation matrix corresponding to the pole of rapol, decpol
    and with zero of new latitude corresponding to ra=ra0
    This matrix need to be np.dot'ed with the input vector to get
    forward transform
    """
    tmppol = np.array(torect(rapol, decpol))  # pole axis
    tmpvec1 = np.array(torect(ra0, 0))  # x axis
    tmpvec1 = np.array(tmpvec1)

    tmpvec1[2] = (-tmppol[0] * tmpvec1[0] - tmppol[1] * tmpvec1[1]) / tmppol[2]
    tmpvec1 /= np.sqrt((tmpvec1**2).sum())
    tmpvec2 = np.cross(tmppol, tmpvec1)  # y axis
    M = np.array([tmpvec1, tmpvec2, tmppol])
    return M


def sphere_rotate(ra, dec, rapol, decpol, ra0, revert=False):
    """ rotate ra,dec to a new spherical coordinate system where the pole is
    at rapol,decpol and the zeropoint is at ra=ra0
    revert flag allows to reverse the transformation
    """

    x, y, z = torect(ra, dec)
    M = rotation_matrix(rapol, decpol, ra0)

    if not revert:
        Axx, Axy, Axz = M[0]
        Ayx, Ayy, Ayz = M[1]
        Azx, Azy, Azz = M[2]
    else:
        Axx, Ayx, Azx = M[0]
        Axy, Ayy, Azy = M[1]
        Axz, Ayz, Azz = M[2]
    xnew = x * Axx + y * Axy + z * Axz
    ynew = x * Ayx + y * Ayy + z * Ayz
    znew = x * Azx + y * Azy + z * Azz
    del x, y, z
    tmp = fromrect(xnew, ynew, znew)
    return (tmp[0], tmp[1])


def rotate_pm(ra, dec, pmra, pmdec, rapol, decpol, ra0, revert=False):
    """
    Rotate the proper motion to the sphere_rotate coord system
    that is specified by the pole direction and right ascencion of the (0,0) pt
    I assume all angles are in degrees and proper motions are in mas/yr
    Arguments:
    ra:
    dec:
    pmra:
    pmdec:
    rapol: float
        RA of the pole
    decpol float
        Dec of the pole
    ra0: float
        RA of the (0,0) point of the new coordinate system
    revert: bool
        if true do the inverse transformation. I.e. convert
        phi1,phi2,pmphii1,pmphi2
        to pmra,pmdec
    Returns:
    pmphi1, pmphi2 in the new coordinate system
    """

    ra, dec, pmra, pmdec = [np.atleast_1d(_) for _ in [ra, dec, pmra, pmdec]]
    M = rotation_matrix(rapol, decpol, ra0)
    if revert:
        M = M.T
    # unit vectors
    e_mura = np.array([-sind(ra), cosd(ra), ra * 0])
    e_mudec = np.array(
        [-sind(dec) * cosd(ra), -sind(dec) * sind(ra),
         cosd(dec)])
    # velocity vector in arbitrary units
    V = pmra * e_mura + pmdec * e_mudec
    del e_mura, e_mudec
    # apply rotation to velocity
    V1 = M @ V
    del V
    X = np.array([cosd(ra) * cosd(dec), sind(ra) * cosd(dec), sind(dec)])
    # apply rotation to position
    X1 = M @ X
    del X
    # rotated coordinates in radians
    lon = np.arctan2(X1[1, :], X1[0, :])
    lat = np.arctan2(X1[2, :], np.sqrt(X1[0, :]**2 + X1[1, :]**2))
    del X1
    # unit vectors in rotated coordinates
    e_mura = np.array([-np.sin(lon), np.cos(lon), lon * 0])
    e_mudec = np.array(
        [-np.sin(lat) * np.cos(lon), -np.sin(lat) * np.sin(lon),
         np.cos(lat)])
    del lon, lat
    return np.sum(e_mura * V1, axis=0), np.sum(e_mudec * V1, axis=0)

def correct_pm(ra, dec, pmra, pmdec, dist):
    """Corrects the proper motion for the speed of the Sun
    Arguments:
        ra - RA in deg
        dec -- Declination in deg
        pmra -- pm in RA in mas/yr (with cosine term)
        pmdec -- pm in declination in mas/yr
        dist -- distance in kpc
    Returns:
        (pmra,pmdec) the tuple with the proper motions corrected for the
        Sun's motion
    """
    C = acoo.ICRS(ra=ra * auni.deg,
                  dec=dec * auni.deg,
                  radial_velocity=0 * kms,
                  distance=dist * auni.kpc,
                  pm_ra_cosdec=pmra * masyr,
                  pm_dec=pmdec * masyr)
    frame = acoo.Galactocentric(**GCPARAMS)
    Cg = C.transform_to(frame)
    Cg1 = acoo.Galactocentric(x=Cg.x,
                              y=Cg.y,
                              z=Cg.z,
                              v_x=Cg.v_x * 0,
                              v_y=Cg.v_y * 0,
                              v_z=Cg.v_z * 0,
                              **GCPARAMS)
    C1 = Cg1.transform_to(acoo.ICRS())
    return ((C.pm_ra_cosdec - C1.pm_ra_cosdec).to_value(masyr),
            (C.pm_dec - C1.pm_dec).to_value(masyr))


def get_CMD_interpolator():
    coloff = -0.01
    magoff = -0.1

    iso_dartmouth_g = np.array([
        14.1916, 13.8004, 13.4031, 12.9995, 12.6317, 12.3286, 12.0868, 11.8791,
        11.6837, 11.4824, 11.2802, 11.0934, 10.9082, 10.7581, 10.6532, 10.586,
        10.5413, 10.5083, 10.4813, 10.4547, 10.426, 10.3941, 10.3567, 10.3109,
        10.2414, 9.9928, 9.7144, 9.5754, 9.4487, 9.3196, 9.1814, 9.0267,
        8.8643, 8.706, 8.5551, 8.4145, 8.2813, 8.15, 8.0148, 7.8739, 7.7934,
        7.7147, 7.637, 7.5598, 7.4838, 7.4089, 7.3364, 7.2656, 7.1963, 7.1287,
        7.0627, 6.9984, 6.9357, 6.8747, 6.8153, 6.7572, 6.7005, 6.645, 6.5907,
        6.5374, 6.4851, 6.4338, 6.3837, 6.3347, 6.2869, 6.2399, 6.1939, 6.1487,
        6.1045, 6.0611, 5.9994, 5.9382, 5.8771, 5.8163, 5.7558, 5.6955, 5.6354,
        5.5754, 5.5156, 5.4559, 5.3965, 5.3375, 5.2791, 5.2213, 5.1642, 5.1079,
        5.0524, 4.9977, 4.9436, 4.8902, 4.8372, 4.7848, 4.7329, 4.6817, 4.6309,
        4.5808, 4.5312, 4.482, 4.4333, 4.385, 4.3699, 4.3544, 4.339, 4.3232,
        4.3073, 4.2913, 4.2752, 4.259, 4.2426, 4.2262, 4.2096, 4.1929, 4.176,
        4.159, 4.1419, 4.1246, 4.1073, 4.0899, 4.0723, 4.0547, 4.0369, 4.0192,
        4.0013, 3.9833, 3.9652, 3.9471, 3.9289, 3.9106, 3.8923, 3.8739, 3.8555,
        3.837, 3.8185, 3.7999, 3.7814, 3.7629, 3.7444, 3.726, 3.7077, 3.6895,
        3.6715, 3.6537, 3.6362, 3.6189, 3.6018, 3.5851, 3.5688, 3.5529, 3.5373,
        3.5228, 3.5145, 3.5072, 3.5006, 3.4933, 3.4841, 3.4738, 3.4632, 3.4515,
        3.4384, 3.4234, 3.4063, 3.387, 3.3656, 3.3425, 3.3181, 3.2925, 3.2661,
        3.2391, 3.2118, 3.184, 3.1526, 3.1202, 3.0872, 3.0534, 3.0192, 2.9844,
        2.9492, 2.9137, 2.8777, 2.8416, 2.805, 2.7678, 2.7301, 2.6921, 2.6533,
        2.6143, 2.5749, 2.5341, 2.4911, 2.4469, 2.4014, 2.3549, 2.3073, 2.2589,
        2.2096, 2.1599, 2.1095, 2.0585, 2.0069, 1.9546, 1.9018, 1.8486, 1.7948,
        1.7407, 1.6863, 1.6768, 1.6668, 1.6566, 1.6464, 1.6365, 1.6264, 1.6162,
        1.6064, 1.5965, 1.5863, 1.474, 1.3721, 1.2678, 1.1626, 1.063, 0.9624,
        0.8616, 0.7646, 0.6683, 0.5721, 0.4812, 0.3935, 0.3142, 0.2533, 0.2509,
        0.2242, 0.1251, 0.0243, -0.0714, -0.165, -0.2543, -0.3415, -0.4262,
        -0.5093, -0.5895, -0.6659, -0.7401, -0.8117, -0.8808, -0.9467, -1.0111,
        -1.0726, -1.1319, -1.1889, -1.2441, -1.2969, -1.3477, -1.3951, -1.4404,
        -1.4833, -1.5249, -1.5657, -1.6049, -1.6429, -1.6793, -1.7145, -1.7484,
        -1.7809, -1.8122, -1.8428, -1.8719, -1.8998, -1.9234
    ])

    iso_dartmouth_r = np.array([
        12.5842, 12.2175, 11.8471, 11.4737, 11.1368, 10.8625, 10.6455, 10.4599,
        10.2854, 10.1058, 9.9258, 9.7602, 9.5962, 9.4637, 9.3714, 9.3122,
        9.273, 9.2441, 9.2205, 9.1974, 9.1726, 9.1447, 9.112, 9.0719, 9.0111,
        8.7983, 8.5647, 8.4486, 8.3433, 8.2368, 8.1238, 7.9988, 7.8693, 7.7443,
        7.6265, 7.5175, 7.4151, 7.3147, 7.2116, 7.1038, 7.0416, 6.9804, 6.9199,
        6.86, 6.8008, 6.7426, 6.6853, 6.629, 6.5737, 6.5194, 6.4663, 6.4143,
        6.3635, 6.3138, 6.2652, 6.2177, 6.1712, 6.1257, 6.081, 6.0372, 5.9942,
        5.952, 5.9105, 5.8698, 5.8298, 5.7905, 5.7517, 5.7137, 5.6763, 5.6395,
        5.587, 5.5345, 5.4821, 5.4298, 5.3774, 5.3252, 5.273, 5.2211, 5.1692,
        5.1172, 5.0653, 5.0136, 4.9622, 4.9112, 4.8607, 4.8107, 4.761, 4.7117,
        4.663, 4.6148, 4.5673, 4.5202, 4.4734, 4.4268, 4.3804, 4.3341, 4.2879,
        4.2418, 4.196, 4.1503, 4.1361, 4.1216, 4.1071, 4.0924, 4.0777, 4.0628,
        4.0478, 4.0327, 4.0174, 4.0018, 3.986, 3.97, 3.9538, 3.9374, 3.9208,
        3.904, 3.8871, 3.87, 3.8528, 3.8354, 3.8178, 3.8001, 3.782, 3.7638,
        3.7454, 3.7268, 3.708, 3.689, 3.6697, 3.6502, 3.6306, 3.6108, 3.5909,
        3.5707, 3.5504, 3.5299, 3.5093, 3.4885, 3.4675, 3.4464, 3.425, 3.4035,
        3.3819, 3.3601, 3.3382, 3.3163, 3.2942, 3.2721, 3.2501, 3.2281, 3.2084,
        3.1879, 3.1675, 3.1473, 3.1269, 3.1062, 3.0853, 3.0637, 3.0411, 3.0169,
        2.9913, 2.9644, 2.9363, 2.9075, 2.8779, 2.8478, 2.8174, 2.7866, 2.7555,
        2.7242, 2.6892, 2.6536, 2.6175, 2.5809, 2.5442, 2.5071, 2.4696, 2.432,
        2.3942, 2.3563, 2.318, 2.2792, 2.2401, 2.2004, 2.1602, 2.1197, 2.079,
        2.0367, 1.9924, 1.9466, 1.8995, 1.8514, 1.8021, 1.752, 1.7009, 1.6492,
        1.5967, 1.5434, 1.4894, 1.4348, 1.3796, 1.3238, 1.2675, 1.2106, 1.1533,
        1.1429, 1.1325, 1.1219, 1.1113, 1.1007, 1.0901, 1.0794, 1.0687, 1.058,
        1.0472, 0.9282, 0.8196, 0.7079, 0.5949, 0.487, 0.3786, 0.2698, 0.165,
        0.0607, -0.0441, -0.1435, -0.2399, -0.3271, -0.394, -0.3956, -0.4237,
        -0.5341, -0.6472, -0.7552, -0.8614, -0.9633, -1.0632, -1.1607, -1.2567,
        -1.3499, -1.4398, -1.5283, -1.6146, -1.6986, -1.7795, -1.8595, -1.9366,
        -2.0117, -2.0847, -2.1561, -2.2251, -2.2923, -2.357, -2.42, -2.4803,
        -2.5395, -2.5976, -2.6541, -2.7092, -2.7626, -2.8146, -2.8653, -2.9145,
        -2.9622, -3.0091, -3.0542, -3.0979, -3.1371
    ])

    CMD_II = scipy.interpolate.UnivariateSpline(
        iso_dartmouth_r[::-1] + magoff,
        (iso_dartmouth_g - iso_dartmouth_r - coloff)[::-1],
        s=0)
    return CMD_II


def read_data():
    # reading data
    from idlsave import idlsave
    import sqlutilpy

    from_cache = True
    if not from_cache:
        mind = 80
        maxd = 100
        Q3 = f'''
    with x as materialized
    (select ebv, ra,dec, flux_g, flux_r, flux_z,flux_w1,flux_ivar_w1,
    type='PSF' as startyp,
    maskbits, anymask_g, anymask_r, anymask_z,
    allmask_g,allmask_r,allmask_z
     from decals_dr9.main as d where dec>-20 and
     q3c_dist(ra,dec,34.5987, 29.7331) between {mind} and {maxd} and flux_r>1
    and ((release=9011 and dec>=32.375) or ( release!=9011 and dec<32.375))
    )
    select * from x left join
    lateral(select q3c_dist(x.ra,x.dec,g.ra,g.dec)*3600 as gdist,
    ruwe, pmra, pmdec, pmra_error,pmdec_error,
    parallax, parallax_error, phot_g_mean_mag, phot_bp_mean_mag,
    phot_rp_mean_mag, source_id,
    pseudocolour,ecl_lat,astrometric_params_solved,nu_eff_used_in_astrometry
    from gaia_dr3.gaia_source
    as g where q3c_join(x.ra,x.dec,g.ra,g.dec,1./3600)
    order by q3c_dist(x.ra,x.dec,g.ra,g.dec) asc limit 1)  as aa  on true ;'''
        print('decals')
        X = sqlutilpy.get(Q3, asDict=True)
        print(len(X['ra']))
        idlsave.save('decals_only_gd1.psav', 'D', X)
    else:
        X, = idlsave.restore('decals_only_gd1.psav', 'D')
    return X


def PM1TRACK(x):
    return scipy.interpolate.UnivariateSpline(
        [-90, -70, -50, -30, -15, 0, 20],
        [-5.5, -6.3, -8, -6.5, -5.7, -5, -3.5])(x)


def PM2TRACK(x):
    return scipy.interpolate.UnivariateSpline([-90, -60, -45, -30, 0, 20],
                                              [-.7, -.7, -0.35, -0.1, 0.2, .5],
                                              s=0)(x)


def pm12_sel_func(gd1fi1, pmfi1, pmfi2, pm_err, pad, mult):
    return np.sqrt((pmfi2 - PM2TRACK(gd1fi1))**2 +
                   (pmfi1 - PM1TRACK(gd1fi1))**2) < pad + mult * pm_err


def DIST(phi1):
    dm = 18.82 + ((phi1 + 48) / 57)**2 - 4.45  # similar to Koposov2010 paper
    return 10**(dm / 5. - 2)


# select to be with X plx errors
def plx_sel_func(gd1fi1, D, mult):
    # extra plx systematic error padding
    plx_sys = 0.05
    subset = np.in1d(D['astrometric_params_solved'], [31, 95])
    plx_zpt_tmp = gaia_zpt.zpt.get_zpt(D['phot_g_mean_mag'][subset],
                                       D['nu_eff_used_in_astrometry'][subset],
                                       D['pseudocolour'][subset],
                                       D['ecl_lat'][subset],
                                       D['astrometric_params_solved'][subset],
                                       _warnings=False)
    plx_zpt = np.zeros(len(D['ra']))
    plx_zpt_tmp[~np.isfinite(plx_zpt_tmp)] = 0
    plx_zpt[subset] = plx_zpt_tmp
    plx = D['parallax'] - plx_zpt
    dplx = 1 / DIST(gd1fi1) - plx

    return np.abs(dplx) < plx_sys + mult * D['parallax_error']


def TRACK(x):
    return scipy.interpolate.CubicSpline(
        [-90, -70, -50, -40, -20, 0, 20],
        [-3, -1.5, -.2, -0., -.0, -1.2, -3])(x)


if __name__ == '__main__':

    D = read_data()

    rapol, decpol, ra_ref = 34.5987, 29.7331, 200

    gd1fi1, gd1fi2 = sphere_rotate(D['ra'], D['dec'], rapol, decpol, ra_ref)

    xpmra, xpmdec = correct_pm(D['ra'], D['dec'], D['pmra'], D['pmdec'],
                               DIST(gd1fi1))
    # reflex correction
    pmfi1, pmfi2 = rotate_pm(D['ra'], D['dec'], xpmra, xpmdec, rapol, decpol,
                             ra_ref)

    pm_err = np.sqrt(0.5 * (D['pmra_error']**2 + D['pmdec_error']**2))

    ext_coeff = dict(g=3.237, r=2.176, z=1.217)
    eg, er, ez = [ext_coeff[_] * D['ebv'] for _ in 'grz']
    ext = {}
    ext['g'] = eg
    ext['r'] = er
    ext['z'] = ez

    g, r, z = [22.5 - 2.5 * np.log10(D['flux_' + _]) - ext[_] for _ in 'grz']

    dfi2 = gd1fi2 - TRACK(gd1fi1)

    CMD_II = get_CMD_interpolator()

    delta_cmd = g - r - CMD_II(r - 5 * np.log10(DIST(gd1fi1) * 1e3) + 5)

    # Selection

    faint_limit = 21
    bright_limit = 16
    # Within the field
    field_sel = betw(dfi2, -10, 10) & betw(gd1fi1, -90, 21)
    # Gaia selection
    pm_pad = 2  # mas/yr padding in pm selection
    gaia_astrom_sel = pm12_sel_func(gd1fi1, pmfi1, pmfi2, pm_err, 
                                    pm_pad, 2.5) & plx_sel_func(
                                        gd1fi1, D, 2.5) & (r > bright_limit)
    print('objects in the field', field_sel.sum())
    print('objects in the field with correct astrometry',
          (gaia_astrom_sel & field_sel).sum())
    bright_cmd_sel = betw(delta_cmd, -.2, .2)  # bigger padding

    stellar_locus_blue_sel = ((betw(r - z - (-.17 + .67 * (g - r)), -0.2, 0.2)
                               & ((g - r) <= 1.1)))
    stellar_locus_red_sel = (((g - r > 1.1)
                              & betw(g - r - (1.05 + .25 * (r - z)), -.2, .2)))
    stellar_locus_sel = stellar_locus_blue_sel | stellar_locus_red_sel

    print('objects in the field with correct astromety and correct cmd',
          (field_sel & gaia_astrom_sel & bright_cmd_sel).sum())

    # no gaia astrometry

    cmd_win = 0.1 + 10**(-2 + (r - 20) / 2.5)
    # 0.11 at r=20 and 0.2 at r=22.5
    faint_sel = ((~np.isfinite(D['pmra'])) & betw(r, 20, faint_limit)
                 & betw(np.abs(delta_cmd), 0,
                        cmd_win)) & D['startyp'] & stellar_locus_sel

    print((faint_sel & field_sel).sum())
    # PSF type + blue in colour and not previously selected
    filler_sel = betw(g - r, -.3, 1.2) & betw(
        r, 19, faint_limit) & D['startyp'] & (~faint_sel) & (
            ~gaia_astrom_sel) & stellar_locus_sel
    filler_red_sel = betw(g - r, 1.2, 2.2) & betw(
        r, 19, faint_limit) & D['startyp'] & (~faint_sel) & (
            ~gaia_astrom_sel) & stellar_locus_sel

    print((filler_sel & field_sel).sum())
    cat1_sel = bright_cmd_sel & gaia_astrom_sel & field_sel
    cat2_sel = faint_sel & field_sel
    cat3_sel = filler_sel & field_sel

    all_cats = cat1_sel | cat2_sel | cat3_sel

    print('Total', all_cats.sum())
    assert (cat1_sel.astype(int) + cat2_sel.astype(int) + cat3_sel.astype(int)
            ).max() <= 1
    PRIORITY = cat1_sel.astype(int) * 10000 + cat2_sel.astype(
        int) * 1000 + cat3_sel.astype(int) * 10
    # ID: int or str, unique for each target.
    # RA, DEC: double, deg
    # SAMPLE: str. This indicates different subsamples of the program that may be
    # assigned different priorities later on. If you only have one sample, give it
    # the name PRINCIPAL.
    # CHECKER: str. Initials identifying person submitting the targets. 2-3 letters
    # PMRA, PMDEC, REF_EPOCH: doubles, I think mas/year, mas/year, Julian year.
    # Set to 0, 0, 2015.5 if not relevant.
    SAMPLE = np.where(
        cat1_sel,
        'BRIGHT_PM',
        np.where(
            cat2_sel,
            'FAINT_NOPM',
            'FILLER'  # np.where(cat3_sel, 'FILLER', 'RED_FILLER')
        ))

    xid = healpy.ang2pix(1 << 29, D['ra'], D['dec'], nest=True,
                         lonlat=True)[all_cats]
    assert (len(np.unique(xid)) == len(xid))
    tabname = 'gd1whole-240109test.fits'
    atpy.Table(
        {
            'ID': xid,
            'RA': D['ra'][all_cats],
            'DEC': D['dec'][all_cats],
            'SAMPLE': SAMPLE[all_cats],

            # 'PRIORITY': PRIORITY[all_cats],
            'FLUX_G': D['flux_g'][all_cats],
            'FLUX_R': D['flux_r'][all_cats],
            'PMRA': D['pmra'][all_cats],
            'PMDEC': D['pmdec'][all_cats],
            'REF_EPOCH': [2016.] * all_cats.sum(),
            'CHECKER': ['SK'] * all_cats.sum()
        },
        meta={
            'CREATED': str(datetime.datetime.now()),
            'CWD': os.getcwd()
        }).write(tabname, overwrite=True)

[geordie666]

@segasai: Thanks! I'll take this code and run with it so I can make it fit some of the desitarget conventions. I'll also add the cross-match between Gaia and DR9 of the Legacy Surveys. I'll also think about how to handle the dependency on gaiadr3_zeropoint.

I'm traveling over the U.S. long weekend. But, at some point next week I'll put together a PR based on this code and send it to you for review.

[segasai]

Great ! Thank you @geordie666 !

Just a small remark, ideally it'd be good if functions rotate_pm, sphere_rotate, correct_pm were not changed (as those are essentially external functions copied here, so if their behaviour was different from what we use elsewhere, it'll be source of pain/confusion). Also if we ever add more streams, likely the basic structure of the selection code will be the same and it'll still rely on those utility functions.

[geordie666]

Added in #814.