![[RSS]](../theme/image/rss.png)
Making Custom Indexes for astrometry.net
When you have an image or a scan of a photographic plate, you usually only have a vague idea of what position the telescope actually was pointed at. Furnishing the image with (more or less) precise information about what pixel correesponds to what sky position is called astrometric calibration. For a while now, arguably the simplest option to do astrometric calibration has been a package called astrometry.net. The eponymous web page has been experiencing… um… operational problems lately, but thanks to the Debian astronomy team, there is a nice package for it in Debian.
However, just running apt install astrometry.net will not give you a working setup. Astrometry.net in addition needs an “index”, files that map star patterns (“quads“, in astrometry.net jargon) to positions. Debian comes with two pre-made sets of indexes at the moment (see apt search astrometry-data): those based on the Tycho 2 catalogue, and those based on 2MASS.
For the index based on Tycho 2, you will find packages astrometry-data-tycho2-10-19, astrometry-data-tycho2-09, astrometry-data-tycho2-08, astrometry-data-tycho2-07[1]. The numbers in there (“scale numbers”) define the size of images the index is good for: 19 means “a major part of the sky”, 10 is “about a degree”, 8 “about half a degree”. Indexes for large images only have a few bright stars and hence are rather compact, which is why 10 though 19 fit into one package, whereas astrometry-data-tycho2-07-littleendian weighs in at 141 MB, and indexes at scale number 0 (suitable for images of a few arcminutes) take dozens of Gigabytes if they are for the whole sky.
So, when you do astrometric calibration, consider the size of your images first and then decide which scale number is sensible for you. It is usually a good idea to try the neighbouring scale numbers, too.
You can then feed these to your calibration routine. If you are running DaCHS, you will probably want to use the AnetHeaderProcessor, where you give the names of the indexes in the sp_indices; you also have to say where to find the indexes, as in:
from gavo import api
class MyObsCalibrator(api.AnetHeaderProcessor):
indexPath = "/usr/share/astrometry"
sp_indices = ["index-tycho2-09*.fits",
"index-tycho2-10*.fits",
"index-tycho2-11*.fits",]
This would be suitable for images that cover about a degree on the sky.
Custom Indexes for Targeted Observations
The Tycho catalogue starts becoming severely incomplete below mV ≈ 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.
| [1] | You will also find that each of these exist in a littleendian and bigendian flavours; ignore these, your machine will pick what it needs when you install the packages without tags. |
