Making Custom Indexes for astrometry.net

Custom Indexes for Targeted Observations

The Tycho catalogue starts becoming severely incomplete below m_V ≈ 11, and since astrometry.net needs a few stars on an image to be able to calibrate it, you cannot use it to calibrate images smaller than a few tens of arcminutes (depending on where you look, of course). If you have smaller images, there are the 2MASS-based indexes; but the bluer your images are, the worse 2MASS as an infrared survey will do, and in addition, having the giant indexes is a big waste of storage and compute resources when you know your images are on a rather small part of the sky.

In such a situation, you will save a lot of CPU and possibly even improve your astrometry if you create a custom index for your specific data. For instance, assume you have images sized about 10 arcminutes, and the observation programme covers a reasonably small set of objects (as long as it's of order a few hundred, a custom index certainly will be a good deal). You could then make your index based on Gaia positions and photometry like this:

"""
Create an index for astrometry.net and a few small fields based on Gaia.

Be sure to adapt this for your use case; for instance, if what your are
calibrating will be from only a part of the sky, pick specific healpixes
(perhaps on a different level; below, we're using level 5).  Also consider
changing the target epoch, the photometry, or the magnitude limit.

This script takes the sample positions from a text file; have
space-separated pairs of ra and dec in targets.txt.
"""

import os
import subprocess

from astropy.table import Table
import pyvo

# 0 is for images of about two arcminutes, 10 for about degree, 12 for two
# degrees, etc.
SIZE_PRESET = 1

# The typical radius of your images in degrees (this is the size of our cone
# searches, so cut some slack); this needs to be changed in unison with
# SIZE_PRESET
IMAGE_RADIUS = 1/10.


def get_target_table():
    """must return an astropy table with columns ra and dec in degrees.

    (of course, if you have your data in a proper format with actual metadata,
    you don't need any of the ugly magic).
    """
    targets = Table.read("targets.txt", format="ascii")
    targets["col1"].name, targets["col2"].name = "ra", "dec"
    targets["ra"].meta = {"ucd": "pos.eq.ra;meta.main"}
    targets["dec"].meta = {"ucd": "pos.eq.dec;meta.main"}
    return targets


def main():
    tap_service = pyvo.dal.TAPService("http://dc.g-vo.org/tap")
    res = tap_service.run_async(f"""
        SELECT g.ra as RA, g.dec as DEC, phot_g_mean_mag as MAG
        FROM gaia.dr3lite AS g
        JOIN TAP_UPLOAD.t1 as mine
            ON DISTANCE(mine.ra, mine.dec, g.ra, g.dec)<{IMAGE_RADIUS}""",
      uploads={"t1": get_target_table()})

    cat_file = "basic-cat.fits"
    res.to_table().write(cat_file, format="fits", overwrite=True)

    try:
        subprocess.run(["build-astrometry-index", "-i", cat_file,
            "-o", f"./index-custom-{SIZE_PRESET:02d}.fits",
            "-P", str(SIZE_PRESET), "-S", "MAG"])
    finally:
        os.unlink(cat_file)


if __name__=="__main__":
    main()

This writes a single file, index-custom-01.fits (in this case).

If you read your positions from something else than the simple ASCII file I'm assuming here: Be sure to annotate the columns containing RA and Dec with the proper UCDs as shown here. That makes DaCHS (and perhaps other TAP services, too) create the right hints for the database, speeding up things tremendously.

You can of course change the ADQL query; it might, for instance, help to replace the G magnitudes with RP or BP ones, or you could use a different catalogue than Gaia. Just make sure the FITS table that is written to basic-cat.fits has exactly the columns RA, DEC, and MAG.

In DaCHS, I tend to keep scripts like the one above in a subdirectory of the resdir called custom-index, and then in the calibration script I write:

from gavo import api

RD = api.getRD("myres/q")

class MyObsCalibrator(api.AnetHeaderProcessor):
  indexPath = RD.resdir
  sp_indices = ["custom-index/index-custom-01.fits"]

Custom Indexes for Ancient Observations

On the other hand, if you have oldish images not going terribly deep, you may want to tailor an index for about the epoch the images were taken at. Many bright stars have a proper motion large enough to matter over a century, and so doing epoch propagation (in this case with the ivo_epoch_prop user defined function, which is not available everywhere) is probably a good idea. The following script computes three full-sky indexes with quads around the desired size; note how you can set the limiting magnitude and the size preset:

"""
Create a full-sky index for bright stars and astrometry.net based on Gaia.

This only works for rather bright stars because the Gaia service will refuse
to server more than ~1e7 objects.

Make sure to choose SIZE_PRESET to your use case (19 means 30 deg,
10 about a degree, two preset steps are about a factor two in scale).
"""

import os
import subprocess

import pyvo

# see the module docstring
SIZE_PRESET = 12

# ignore stars fainter than this; you can't go below 14 all-sky with Gaia
# and the GAVO DC server
MAX_MAG = 12

# Epoch to transform the stars to
TARGET_EPOCH = 1910


def main():
    tap_service = pyvo.dal.TAPService("http://dc.g-vo.org/tap")
    res = tap_service.run_async(f"""
        SELECT pos[1] as RA, pos[2] as DEC, mag as MAG
        FROM (
            SELECT phot_bp_mean_mag AS mag,
                ivo_epoch_prop(ra, dec, parallax,
                    pmra, pmdec, radial_velocity, 2016, {TARGET_EPOCH}) as pos
            FROM gaia.dr3lite
          WHERE phot_bp_mean_mag<{MAX_MAG}) AS q""")

    cat_file = "current.fits"
    res.to_table().write(cat_file, format="fits", overwrite=True)

    try:
        for size_preset in range(SIZE_PRESET-1, SIZE_PRESET+2):
            subprocess.run(["build-astrometry-index", "-i", cat_file,
                "-o", f"./index-custom-{size_preset:02d}.fits",
                "-P", str(size_preset), "-S", "MAG"])
    finally:
        os.unlink(cat_file)


if __name__=="__main__":
    main()

With this and my custom-index directory, your DaCHS header processor could say:

from gavo import api

RD = api.getRD("myres/q")

class MyObsCalibrator(api.AnetHeaderProcessor):
  indexPath = RD.resdir
  sp_indices = ["custom-index/index-custom-*.fits"]

Custom Indexes: Full-sky and Deep

I have covered the cases “deep and spotty” and “shallow and full-sky“. The case “deep and full-sky“ is a bit more involved because it still lies in the realm of big data, which always requires extra tricks. In this case, that would be retrieving the basic catalogue in parts – for instance, by HEALPix – and at the same time splitting the index up between HEALPixes, too. This does not require great magic, but it does require a bit of non-trivial bookkeeping, and hence I will only write about it if someone actually needs it – if that's you, please write in.