Astrometric calibration#
This page presents a way get an accurate WCS-transformation for the science frames.
We make use of the local astrometry.net (see Requirements and installation), calling it from within python via subprocess. The executable of astrometry.net that we will use is solve-field. For each image, we will call solve-field twice, with different settings. The first iteration is for rough plate-solving (identifying where the image is on the sky), and the second iteration is to fine-tune a distortion model. The page Requirements and installation gives more information and links about astrometry.net. You are of course very welcome to learn more about this rather fascinating tool, but for the purpose of this lab course this notebook could in principle be run as a black box.
The ambition of this notebook is to obtain a WCS that is accurate enough for forced aperture photometry or stacking, i.e., with astrometric residuals significantly smaller than a pixel. This notebook will therefore also create a plot of astrometric residuals of each image to check this.
One note on the algorithmic approach: all our images are taken with the same camera/telescope, however we will constrain the distortion model independently for every image. While this might not seem to be very elegant (in particular given a potentially low number of stars per image), this approach remains very flexible, and works even when the camera got rotated (for example) or when flexion in the telescope affects the distortion.
Note
In case you’re using this tutorial with images from a different telescope setup, simply either adapt (or remove!) the hard-coded pixel scale parameters --scale-low and --scale-high as well as --downsample from the call to solve-field below. These parameters are optional, they are provided here only to speed up the process. Note that if your telescope has a field of view much narrower or wider than about 0.5 deg, you might need different index files (see links on Requirements and installation).
As before, you could copy or write the code shown below in a script, or alternatively directly download this page as a jupyter notebook file.
To run the code, you’ll need the module dataredconfig.py, as explained here.
import dataredconfig
import subprocess
import pathlib
import numpy as np
import astropy
import astropy.visualization
import astropy.table
from astropy import units as u
%matplotlib widget
import matplotlib
from matplotlib import pyplot as plt
import ccdproc
# Locating the image files to run on:
input_dir = dataredconfig.work_dir / "LIGHT_PRERED"
science_files = ccdproc.ImageFileCollection(input_dir, keywords=dataredconfig.ifc_header_keywords)
# You may want to use only some of the files in there, for example:
#object_to_process = "HD92670"
#science_files = science_files.filter(object=object_to_process)
# Overview of science files:
science_files.summary
# You may want to try this instead, for a better display:
#science_files.summary.show_in_notebook()
# Directory in which to write the output files:
astrometry_dir = dataredconfig.work_dir / "ASTROMETRY"
astrometry_dir.mkdir(exist_ok=True)
# Subdirectories for the two runs:
run_1_dir = astrometry_dir / "run_1"
run_1_dir.mkdir(exist_ok=True)
run_2_dir = astrometry_dir # The output of run2 is directly written into the astrometry_dir
run_2_dir.mkdir(exist_ok=True)
# So after this notebook is done, you will want to use the files that are directly inside your astrometry_dir for the next steps.
We now define one function for visualization. As you will see once it has run, the generated plot is relatively self-explanatory.
def visualize_corr_file(corr_file_path, fig_file_path=None, ax=None):
"""Function to vizualize the quality of the astrometric solution
corr_file_path is the path to a .corr file as written by solve-field
Specify either fig_file_path (figure will be written there), or a matplotlib
axes object (will plot on these axes).
"""
# solve-field writes these "corr" tables with all "corresponding" stars it uses, we open them:
corr = astropy.table.Table.read(corr_file_path, format='fits')
print(f"Read {len(corr)} stars from corr file...")
corr["error_x"] = corr["field_x"] - corr["index_x"]
corr["error_y"] = corr["field_y"] - corr["index_y"]
corr["error"] = np.hypot(corr["error_x"], corr["error_y"])
if fig_file_path:
fig, ax = plt.subplots(figsize=(9, 6))
elif not ax:
print("Provide either fig_file_path or ax!")
#ax.plot(corr["field_x"], corr["field_y"], "r.")
q = ax.quiver(corr["field_x"], corr["field_y"], corr["error_x"], corr["error_y"], units='dots', scale=0.02)
qk = ax.quiverkey(q, 0.4, 0.9, 1, "Astrometric error of 1 pixel", labelpos='E', coordinates='figure')
cm = matplotlib.colormaps['RdYlBu_r']
ax.scatter(corr["field_x"], corr["field_y"],
s = 15,
c = corr["error"].value,
cmap=cm
)
ax.set_aspect("equal")
ax.set_xlabel("x [pixel]")
ax.set_ylabel("y [pixel]")
if fig_file_path:
fig.savefig(fig_file_path)
plt.close()
And a function to run solve-field. See the documentation of astrometry.net (http://astrometry.net/doc/readme.html) and the man-page of solve-field (https://manpages.debian.org/bullseye/astrometry.net/solve-field.1.en.html) to learn more about the options we use below.
def run_solvefield(input_fits_path, dir_path, settings=1, input_wcs_path=None):
"""Wrapper to run solve-field in a subprocess, with some hard-coded settings.
The stdout and stderr are written to text files, as "logs" to check in case of problems.
Use settings=2 and provide the wcs of run 1 when running for the second time, to fine-tune the SIP.
A log file gets written along the other output files from solve-field
"""
cmd_1 = ["solve-field",
"--overwrite",
"--scale-units", "arcsecperpix",
"--scale-low", "0.60",
"--scale-high", "0.62",
"--downsample", "4",
"--crpix-center",
"--no-verify", # We ignore existing WCS
"--tweak-order", "2",
"--resort",
"--dir", dir_path,
"--no-plots",
"--new-fits", "%s.fits",
input_fits_path
]
cmd_2 = ["solve-field",
"--overwrite",
"--downsample", "4",
"--crpix-center",
"--tweak-order", "3", # SIP order 3!
"--resort",
"--dir", dir_path,
"--verify", input_wcs_path, # We verify the WCS from run 1
"--no-verify-uniformize",
"--objs", "1000000", # We want to use all available stars
"--no-plots",
"--new-fits", "%s.fits",
input_fits_path
]
if settings == 1:
cmd = cmd_1
elif settings == 2:
cmd = cmd_2
else:
raise RuntimeError("settings must be 1 or 2")
# Creating a log file:
log_filepath = pathlib.Path(dir_path) / pathlib.Path(input_fits_path).with_suffix(".log").name
log_file = open(log_filepath, "w")
# And running solve-field:
res = subprocess.run(cmd, text=True, stdout=log_file, stderr=subprocess.STDOUT)
log_file.close()
print(f"Log file: {log_filepath}")
return(res.returncode)
# We'll ignore some astropy warnings that get raised as the corr files written by solve-field have non-compliant units:
import warnings
warnings.simplefilter('ignore', category=astropy.units.UnitsWarning)
Finally we run over all the input files. In case of a failure (in practice: lack of stars due to clouds etc) the image will simply get ignored.
# Looping over the input files:
for (i, science_file) in enumerate(science_files.files_filtered(include_path=True)):
n_files = len(science_files.summary)
filename_stem = pathlib.Path(science_file).stem
print(f"=== Processing image {i+1}/{n_files}: {filename_stem} ===")
# Run 1:
print(f"Run 1")
try:
run_solvefield(science_file, run_1_dir, settings=1)
except:
print("It failed!")
else:
# Making a checkplot:
corr_file_path = pathlib.Path(run_1_dir) / pathlib.Path(science_file).with_suffix(".corr").name
fig_file_path = pathlib.Path(run_1_dir) / pathlib.Path(science_file).with_suffix(".checkplot.png").name
if corr_file_path.exists():
visualize_corr_file(corr_file_path, fig_file_path=fig_file_path)
print(f"Checkplot: {fig_file_path}")
# Run 2:
print(f"Run 2")
run_2_input_path = pathlib.Path(run_1_dir) / pathlib.Path(science_file).name
run_2_input_wcs_path = pathlib.Path(run_1_dir) / pathlib.Path(science_file).with_suffix(".wcs").name
try:
run_solvefield(run_2_input_path, run_2_dir, settings=2, input_wcs_path=run_2_input_wcs_path)
except:
print("It failed!")
else:
# Making a checkplot:
corr_file_path = pathlib.Path(run_2_dir) / pathlib.Path(science_file).with_suffix(".corr").name
fig_file_path = pathlib.Path(run_2_dir) / pathlib.Path(science_file).with_suffix(".checkplot.png").name
if corr_file_path.exists():
visualize_corr_file(corr_file_path, fig_file_path=fig_file_path)
print(f"Checkplot: {fig_file_path}")