Details on astrometric calibration

%autoawait off
%matplotlib inline
%run 01_notebook_configuration.py

The header of raw FITS files does not usually contain a valid WCS. The astrometric calibration is implemented automatically during the  ReducedFit.build_file() step, and can be done separately through ReducedFit.astrometric_calibration().

from iop4lib.db import ReducedFit
redf = ReducedFit.by_fileloc("OSN-T090/2022-09-18/BLLac-0001R.fit")
redf.build_file()

2024-09-24 16:02 - pid 4128 - [instrument.py:480] - DEBUG - <ReducedFit 145 | OSN-T090/2022-09-18/BLLac-0001R.fit>: building file
2024-09-24 16:02 - pid 4128 - [instrument.py:386] - DEBUG - <ReducedFit 145 | OSN-T090/2022-09-18/BLLac-0001R.fit>: applying masters
2024-09-24 16:02 - pid 4128 - [instrument.py:395] - WARNING - <ReducedFit 145 | OSN-T090/2022-09-18/BLLac-0001R.fit>: no masterdark found, assuming dark current = 0, is this a CCD camera and it's cold?
2024-09-24 16:02 - pid 4128 - [instrument.py:486] - DEBUG - <ReducedFit 145 | OSN-T090/2022-09-18/BLLac-0001R.fit>: performing astrometric calibration
2024-09-24 16:02 - pid 4128 - [astrometry.py:179] - DEBUG - <ReducedFit 145 | OSN-T090/2022-09-18/BLLac-0001R.fit>: 27 different combinations of parameters to try.
2024-09-24 16:02 - pid 4128 - [astrometry.py:182] - DEBUG - <ReducedFit 145 | OSN-T090/2022-09-18/BLLac-0001R.fit>: attempt 1 / 27, ({'keep_n_seg': 150, 'border_margin_px': 5, 'output_logodds_threshold': 14, 'n_rms_seg': 6.0, 'bkg_filter_size': 11, 'bkg_box_size': 16, 'seg_fwhm': 1.0, 'seg_kernel_size': None, 'npixels': 32, 'allsky': False, 'position_hint': PositionHint(ra_deg=330.7291666666667, dec_deg=42.26611111111111, radius_deg=0.32999999999999996)}) ...
2024-09-24 16:02 - pid 4128 - [astrometry.py:275] - DEBUG - <ReducedFit 145 | OSN-T090/2022-09-18/BLLac-0001R.fit>: len(pos_seg)=187
2024-09-24 16:02 - pid 4128 - [astrometry.py:278] - DEBUG - <ReducedFit 145 | OSN-T090/2022-09-18/BLLac-0001R.fit>: Keeping only 150 brightest segments.
2024-09-24 16:02 - pid 4128 - [astrometry.py:282] - DEBUG - <ReducedFit 145 | OSN-T090/2022-09-18/BLLac-0001R.fit>: Removing segments within 5 px from border.
2024-09-24 16:02 - pid 4128 - [astrometry.py:312] - DEBUG - <ReducedFit 145 | OSN-T090/2022-09-18/BLLac-0001R.fit>: trying to solve astrometry with Seg Pos (n=147) (output_logodds_threshold=14).
2024-09-24 16:02 - pid 4160 - [__init__.py:212] - INFO - loaded 240 index files
2024-09-24 16:02 - pid 4160 - [__init__.py:276] - INFO - solve 1: start
2024-09-24 16:02 - pid 4160 - [__init__.py:276] - INFO - solve 1: slice=[0, 25[ (1 / 5), index="5200/index-5200-14.fits" (1 / 5)
2024-09-24 16:02 - pid 4160 - [__init__.py:276] - INFO - solve 1: logodds=380.916, matches=141, conflicts=0, distractors=2, ra=330.695, dec=42.2697, scale=0.773671, index="5200/index-5200-14.fits"
2024-09-24 16:02 - pid 4128 - [astrometry.py:327] - DEBUG - <ReducedFit 145 | OSN-T090/2022-09-18/BLLac-0001R.fit>: Seg Pos (n=147) worked.
2024-09-24 16:02 - pid 4128 - [astrometry.py:328] - DEBUG - <ReducedFit 145 | OSN-T090/2022-09-18/BLLac-0001R.fit>: bm.index_path=PosixPath('/mnt/astrometry_cache/5200/index-5200-14.fits')
2024-09-24 16:02 - pid 4128 - [astrometry.py:329] - DEBUG - <ReducedFit 145 | OSN-T090/2022-09-18/BLLac-0001R.fit>: bm.center_ra_deg=330.6954314398751
2024-09-24 16:02 - pid 4128 - [astrometry.py:330] - DEBUG - <ReducedFit 145 | OSN-T090/2022-09-18/BLLac-0001R.fit>: bm.center_dec_deg=42.26972920042598
2024-09-24 16:02 - pid 4128 - [astrometry.py:331] - DEBUG - <ReducedFit 145 | OSN-T090/2022-09-18/BLLac-0001R.fit>: bm.scale_arcsec_per_pixel=0.7736712092056637
2024-09-24 16:02 - pid 4128 - [astrometry.py:332] - DEBUG - <ReducedFit 145 | OSN-T090/2022-09-18/BLLac-0001R.fit>: bm.logodds=380.9163513183594
2024-09-24 16:02 - pid 4128 - [astrometry.py:192] - DEBUG - <ReducedFit 145 | OSN-T090/2022-09-18/BLLac-0001R.fit>: WCS built with attempt 1 / 27 ({'keep_n_seg': 150, 'border_margin_px': 5, 'output_logodds_threshold': 14, 'n_rms_seg': 6.0, 'bkg_filter_size': 11, 'bkg_box_size': 16, 'seg_fwhm': 1.0, 'seg_kernel_size': None, 'npixels': 32, 'allsky': False, 'position_hint': PositionHint(ra_deg=330.7291666666667, dec_deg=42.26611111111111, radius_deg=0.32999999999999996)}).
2024-09-24 16:02 - pid 4128 - [astrometry.py:204] - DEBUG - <ReducedFit 145 | OSN-T090/2022-09-18/BLLac-0001R.fit>: building summary images.
2024-09-24 16:02 - pid 4128 - [plotting.py:476] - DEBUG - <ReducedFit 145 | OSN-T090/2022-09-18/BLLac-0001R.fit>: plotting astrometry summary image of background substraction results
2024-09-24 16:02 - pid 4128 - [plotting.py:488] - DEBUG - <ReducedFit 145 | OSN-T090/2022-09-18/BLLac-0001R.fit>: plotting astrometry summary image of segmentation results
2024-09-24 16:02 - pid 4128 - [plotting.py:582] - DEBUG - <ReducedFit 145 | OSN-T090/2022-09-18/BLLac-0001R.fit>: plotting astrometry summary image of astrometry results
2024-09-24 16:02 - pid 4128 - [plotting.py:214] - DEBUG - <ReducedFit 145 | OSN-T090/2022-09-18/BLLac-0001R.fit>: found 5 catalog sources in field: [AstroSource.objects.get(name='2200+420'), AstroSource.objects.get(name='BL Lacertae C'), AstroSource.objects.get(name='BL Lacertae H'), AstroSource.objects.get(name='BL Lacertae test 2'), AstroSource.objects.get(name='BL Lacertae test 4')]
2024-09-24 16:02 - pid 4128 - [instrument.py:425] - DEBUG - <ReducedFit 145 | OSN-T090/2022-09-18/BLLac-0001R.fit>: saving WCSs to FITS header.
2024-09-24 16:02 - pid 4128 - [instrument.py:497] - DEBUG - <ReducedFit 145 | OSN-T090/2022-09-18/BLLac-0001R.fit>: astrometric calibration was successful.
2024-09-24 16:02 - pid 4128 - [instrument.py:500] - DEBUG - <ReducedFit 145 | OSN-T090/2022-09-18/BLLac-0001R.fit>: searching for sources in field...
2024-09-24 16:02 - pid 4128 - [instrument.py:503] - DEBUG - <ReducedFit 145 | OSN-T090/2022-09-18/BLLac-0001R.fit>: found 5 sources in field.

For simple photometry images such as those from the Andor cameras, tries to detect sources in the image with different sets of parameters and feeds them to a local astrometry.net solver. It tries with several sets until it successfully calibrates the image. When done through build_file(), if calibration is not achieved, it will give the ERROR_ASTROMETRY flag, otherwise it will give it the BUILT_REDUCED flag.

redf

ReducedFit(id=145):
- telescope: OSN-T090
- night: 2022-09-18
- filename: BLLac-0001R.fit
- instrument: AndorT90
- imgtype: LIGHT
- size: 1024x1024
- obsmode: PHOTOMETRY
- band: R
- exptime: 60.0
- flags: DOWNLOADED,CLASSIFIED,BUILT_REDUCED

redf.flag_labels

['DOWNLOADED', 'CLASSIFIED', 'BUILT_REDUCED']

In the reduction log

2024-02-21 11:51 - pid 75152 - [plotting.py:476] - DEBUG - <ReducedFit 212 | OSN-T090/2022-09-18/BLLac-0001R.fit>: plotting astrometry summary image of background substraction results
2024-02-21 11:51 - pid 75152 - [plotting.py:488] - DEBUG - <ReducedFit 212 | OSN-T090/2022-09-18/BLLac-0001R.fit>: plotting astrometry summary image of segmentation results
2024-02-21 11:51 - pid 75152 - [plotting.py:582] - DEBUG - <ReducedFit 212 | OSN-T090/2022-09-18/BLLac-0001R.fit>: plotting astrometry summary image of astrometry results

you can see that several summary images were created. They are displayed in the admin page for the reduced fit

For images with ordinary and extraordinary pairs such as those from CAFOS and DIPOL, the reduction procedure first finds the sources and then divides them in ordinaty and extraordinary lists. This is done through the utilities at iop4lib.utils.sourcepairing submodule.

To find the pairs, the algorithm builds a distribution of the distance between all sources in the image, and selectes the most frequent distance with some resolution. This process is not free from error. It can be improved by providing some bounds for this expected distance. For example, in a CAFOS image such as

from iop4lib.db import RawFit
rf = RawFit.by_fileloc("CAHA-T220/2022-09-18/caf-20220918-23:01:33-sci-agui.fits")

We first create the reduced image and apply the calibration frames

redf = ReducedFit.create(rawfit=rf)
redf.apply_masters()

2024-09-24 16:02 - pid 4128 - [reducedfit.py:121] - DEBUG - DB entry of ReducedFit for CAHA-T220/2022-09-18/caf-20220918-23:01:33-sci-agui.fits already exists, it will be used instead.
2024-09-24 16:02 - pid 4128 - [instrument.py:373] - WARNING - <ReducedFit 151 | CAHA-T220/2022-09-18/caf-20220918-23:01:33-sci-agui.fits>: masterbias in this epoch could not be found, attemptying adjacent epochs.
2024-09-24 16:02 - pid 4128 - [instrument.py:373] - WARNING - <ReducedFit 151 | CAHA-T220/2022-09-18/caf-20220918-23:01:33-sci-agui.fits>: masterflat in this epoch could not be found, attemptying adjacent epochs.
2024-09-24 16:02 - pid 4128 - [instrument.py:386] - DEBUG - <ReducedFit 151 | CAHA-T220/2022-09-18/caf-20220918-23:01:33-sci-agui.fits>: applying masters
2024-09-24 16:02 - pid 4128 - [instrument.py:395] - WARNING - <ReducedFit 151 | CAHA-T220/2022-09-18/caf-20220918-23:01:33-sci-agui.fits>: no masterdark found, assuming dark current = 0, is this a CCD camera and it's cold?

The .has_pairs property tells whether an image is expected to have pairs:

redf.has_pairs

True

Since we will be doing “blind” astrometric calibration with the local astrometry.net solver, we will need (or at least greatly benefit) from hints for the size of the pixel and for pointing of the camera. These can be obtained with

size_hint = redf.get_astrometry_size_hint()
size_hint

SizeHint(lower_arcsec_per_pixel=0.5035, upper_arcsec_per_pixel=0.5565000000000001)

This code to build this size hint is implemented in the corresponding instrument class.

We will also benefit from a hint for the position of the object in the image.

The position is hint is built by the corresponding instrument class by looking at the RA and DEC (or equivalent) keywords in the FITS header. If the FITS header contains the RA and DEC (or equivalent) keywords, the header_hintcoord property will return an SkyCoord object with the position

redf.header_hintcoord

<SkyCoord (ICRS): (ra, dec) in deg
    (330.9182, 42.38482)>

This position should be close to the observed object. The header_hintobject property returns the AstroSource in the catalog that matches the name in the OBJECT (or equivalent) keyword:

redf.header_hintobject

AstroSource(name='2200+420'):
- other_names: BL Lacertae
- ra=330.6804 deg, dec=42.2777 deg
- srctype: blazar
- comment:

References:

• MAPCAT catalog. OPERACAT table blazar_source.

redf.header_hintobject.coord.separation(redf.header_hintcoord)

\[0^\circ12{}^\prime21.1434803{}^{\prime\prime}\]

It can be be seen that the separation is indeed small.

When header_hintcoord is available, get_astrometry_position_hint() will return a hint for the object built from it. Otherwise, it will look at the header_hintobject. If none are available, it will raise an exception.

position_hint = redf.get_astrometry_position_hint()
position_hint

PositionHint(ra_deg=330.9182, dec_deg=42.38482, radius_deg=0.4)

All these hints will be used to accelerate the blind calibration.

Now, let’s detect the stars visible in the image. The whole process will depend on the parameters we select, and we might need to try with several combinations.

bkg_filter_size = 11 # a higher number will produce a more smooth bkg
bkg_box_size = 16 # a higher number will produce a more detailed bkg
# these bkg parameters should be set so the computed bkg does not include light from the point sources, but is detailed enough to show real variations along the image.
seg_kernel_size = None # none will automatically build a kernel of size 2n+1, which should be almost always a good approximation
npixels = 32 # related to the size of the sources in the image
seg_fwhm = 1.0 # the image will be convolved with a gaussian kernel of this width before src detection
n_rms_seg = 1.0 # the threshold of signal over background to count as a detection
# these detection parameters should return most of the real points in the image without too much fake sources
keep_n_seg = 200 # the maximum number of detected sources to use in the calibration
border_margin_px = 20 # remove any pixels less than this distance from the borders

We can use the utility functions in iop4lib.utils.sourcedetection to get an bkg-estimation, which utilize the photutils package (you could do the same, or try different methods following its docs.

Background substraction

from iop4lib.utils.sourcedetection import get_bkg, get_segmentation, get_cat_sources_from_segment_map
bkg = get_bkg(redf.mdata, filter_size=bkg_filter_size, box_size=bkg_box_size)
imgdata_bkg_substracted = redf.mdata - bkg.background

from iop4lib.utils.plotting import plot_preview_background_substraction_1row
import matplotlib.pyplot as plt
import matplotlib as mplt
fig = plt.figure(figsize=(12, 5))
plot_preview_background_substraction_1row(redf, bkg, fig=fig)
plt.show()

Image Segmentation

seg_threshold = n_rms_seg * bkg.background_rms
segment_map, convolved_data = get_segmentation(imgdata_bkg_substracted, fwhm=seg_fwhm, kernel_size=seg_kernel_size, npixels=npixels, threshold=seg_threshold)
seg_cat, pos_seg, tb = get_cat_sources_from_segment_map(segment_map, imgdata_bkg_substracted, convolved_data)

from astropy.visualization import SqrtStretch, LogStretch, AsymmetricPercentileInterval
from astropy.visualization.mpl_normalize import ImageNormalize
fig, ax = plt.subplots()
ax.imshow(imgdata_bkg_substracted, cmap='gray', norm=ImageNormalize(imgdata_bkg_substracted.compressed(), stretch=LogStretch(a=10), interval=AsymmetricPercentileInterval(30, 99)))
seg_cat.plot_kron_apertures(color='r', lw=1.5)
plt.show()

Source pairing

Now we need to distinguish between ordinary and extraordinary sources. To facilitate this, you can use the get_pairs_dxy function, which will try to separate them by looking at the most common distance between sources.

First we find the most common distance between sources (both in x and y)

from iop4lib.utils.sourcepairing import get_pairs_dxy, get_best_pairs
bins = 60
hist_range = (0, 60)
dx_min = dy_min = 0
dx_max = dy_max = 60
dx_eps = dy_eps = 1.0
seg1xy, seg2xy, seg_disp_xy, seg_disp_sign_xy = get_pairs_dxy(pos_seg, bins=bins, hist_range=hist_range, dx_min=dx_min, dx_max=dx_max, dy_min=dy_min, dy_max=dy_max, dx_eps=dx_eps, dy_eps=dy_eps, doplot=True)
seg_disp_sign_xy

array([-36.33263958,  -0.12733392])

Then we select the best pairs according to this distance, and recalculate it

seg1xy_best, seg2xy_best, seg_disp_xy_best, seg_disp_sign_xy_best = get_best_pairs(seg1xy, seg2xy, seg_disp_sign_xy)
seg_disp_sign_xy_best

array([-36.33263958,  -0.12733392])

We can plot the image again and check that indeed the pairs are well separated

from iop4lib.utils.plotting import imshow_w_sources
imshow_w_sources(redf.mdata, pos1=seg1xy_best, pos2=seg2xy_best)

To avoid mistakes, we provided lower and upper bounds for the distance between pairs. This can be specially important if there are few sources in the image. We could also do the same by ourselves, by looking at the distribution of distances:

import itertools
pairs = list(itertools.combinations(pos_seg, 2))
len(pairs)

35511

dx_min = 0
dx_max = 60
dy_min = 0
dy_max = 60
paired = [(p1,p2) for p1,p2 in pairs if ( dx_min < abs( p1[0] - p2[0] ) < dx_max and dy_min < abs( p1[1] - p2[1] ) < dy_max )]
len(paired)

867

fig, ax = plt.subplots()
for i in [0,1]:
    distances = [abs(p1[i]-p2[i]) for p1,p2 in paired]
    plt.hist(distances, bins=60, range=(0,60), alpha=0.7)
plt.xlabel('Distance (px)')
plt.ylabel('Number of pairs')
plt.show()

Finally, we need to invoke the local astrometry.net solver to calibrate the image.

from iop4lib.utils.astrometry import solve_astrometry

solution = solve_astrometry(seg1xy_best, size_hint=size_hint, position_hint=position_hint)
solution

2024-09-24 16:02 - pid 4197 - [__init__.py:212] - INFO - loaded 240 index files
2024-09-24 16:02 - pid 4197 - [__init__.py:276] - INFO - solve 1: start
2024-09-24 16:02 - pid 4197 - [__init__.py:276] - INFO - solve 1: slice=[0, 25[ (1 / 3), index="5200/index-5200-14.fits" (1 / 5)
2024-09-24 16:02 - pid 4197 - [__init__.py:276] - INFO - solve 1: slice=[0, 25[ (1 / 3), index="5200/index-5201-14.fits" (2 / 5)
2024-09-24 16:03 - pid 4197 - [__init__.py:276] - INFO - solve 1: slice=[0, 25[ (1 / 3), index="5200/index-5202-14.fits" (3 / 5)
2024-09-24 16:03 - pid 4197 - [__init__.py:276] - INFO - solve 1: slice=[0, 25[ (1 / 3), index="5200/index-5203-14.fits" (4 / 5)
2024-09-24 16:03 - pid 4197 - [__init__.py:276] - INFO - solve 1: slice=[0, 25[ (1 / 3), index="5200/index-5204-14.fits" (5 / 5)
2024-09-24 16:03 - pid 4197 - [__init__.py:276] - INFO - solve 1: slice=[25, 50[ (2 / 3), index="5200/index-5200-14.fits" (1 / 5)
2024-09-24 16:03 - pid 4197 - [__init__.py:276] - INFO - solve 1: logodds=160.586, matches=70, conflicts=0, distractors=14, ra=330.675, dec=42.2686, scale=0.524084, index="5200/index-5200-14.fits"

Solution(solve_id='1', matches=[Match(logodds=160.58602905273438, center_ra_deg=330.67545493255574, center_dec_deg=42.26859969749991, scale_arcsec_per_pixel=0.5240836722782216, index_path=PosixPath('/mnt/astrometry_cache/5200/index-5200-14.fits'), stars=(Star(ra_deg=330.689884271765, dec_deg=42.3196494803186, metadata={}), Star(ra_deg=330.7250235380827, dec_deg=42.22172265671082, metadata={}), Star(ra_deg=330.68921881915674, dec_deg=42.276504127886966, metadata={}), Star(ra_deg=330.6686605321736, dec_deg=42.244955317312844, metadata={}), Star(ra_deg=330.606496663498, dec_deg=42.29642927956227, metadata={}), Star(ra_deg=330.62824559333745, dec_deg=42.28602926215711, metadata={}), Star(ra_deg=330.6803806918542, dec_deg=42.27777212498205, metadata={}), Star(ra_deg=330.66214717435116, dec_deg=42.31664362093201, metadata={}), Star(ra_deg=330.668073572888, dec_deg=42.316474287620046, metadata={}), Star(ra_deg=330.63916131581885, dec_deg=42.24035505186057, metadata={}), Star(ra_deg=330.6754327942969, dec_deg=42.21822059373048, metadata={}), Star(ra_deg=330.65423699698255, dec_deg=42.2271577292143, metadata={}), Star(ra_deg=330.62742302432804, dec_deg=42.31657787928523, metadata={}), Star(ra_deg=330.74252280208543, dec_deg=42.305447819823215, metadata={}), Star(ra_deg=330.5991795813756, dec_deg=42.28808295893054, metadata={}), Star(ra_deg=330.72365667038076, dec_deg=42.29682475851014, metadata={}), Star(ra_deg=330.6652497702822, dec_deg=42.227384400395465, metadata={}), Star(ra_deg=330.6120815972558, dec_deg=42.29708764300834, metadata={}), Star(ra_deg=330.6921626029343, dec_deg=42.309079199780356, metadata={}), Star(ra_deg=330.64603270424135, dec_deg=42.313439968191226, metadata={}), Star(ra_deg=330.6775600706942, dec_deg=42.21222599485241, metadata={}), Star(ra_deg=330.6312907472159, dec_deg=42.27082138820124, metadata={}), Star(ra_deg=330.73072341459044, dec_deg=42.230007482225105, metadata={}), Star(ra_deg=330.7386034904528, dec_deg=42.2950819139837, metadata={}), Star(ra_deg=330.6359268542213, dec_deg=42.27984343890986, metadata={}), Star(ra_deg=330.6619423260062, dec_deg=42.24417207647909, metadata={}), Star(ra_deg=330.6037081020162, dec_deg=42.27226594872766, metadata={}), Star(ra_deg=330.67444555564566, dec_deg=42.21723403323365, metadata={}), Star(ra_deg=330.6724467550032, dec_deg=42.31685033592088, metadata={}), Star(ra_deg=330.6503069580795, dec_deg=42.28163176890514, metadata={}), Star(ra_deg=330.63197049203325, dec_deg=42.211112334865646, metadata={}), Star(ra_deg=330.6155332086872, dec_deg=42.246940493628095, metadata={}), Star(ra_deg=330.71458310461696, dec_deg=42.25858707453221, metadata={}), Star(ra_deg=330.7387326142312, dec_deg=42.226806068016636, metadata={}), Star(ra_deg=330.74129924189464, dec_deg=42.28111027710141, metadata={}), Star(ra_deg=330.72791207724276, dec_deg=42.23117127346185, metadata={}), Star(ra_deg=330.65100319670586, dec_deg=42.3087516501433, metadata={}), Star(ra_deg=330.70126242336914, dec_deg=42.312856191395724, metadata={}), Star(ra_deg=330.61572174801864, dec_deg=42.300805296476895, metadata={}), Star(ra_deg=330.68279436072703, dec_deg=42.25795333781425, metadata={}), Star(ra_deg=330.6922200979271, dec_deg=42.22810816105225, metadata={}), Star(ra_deg=330.7178450712782, dec_deg=42.278323519120285, metadata={}), Star(ra_deg=330.7527689297072, dec_deg=42.24973411878969, metadata={}), Star(ra_deg=330.74806130640457, dec_deg=42.266961613056495, metadata={}), Star(ra_deg=330.7017402283431, dec_deg=42.247357107459706, metadata={}), Star(ra_deg=330.63009759101675, dec_deg=42.261163326726255, metadata={}), Star(ra_deg=330.6626570982849, dec_deg=42.31590787406957, metadata={}), Star(ra_deg=330.6911891847441, dec_deg=42.22778913892858, metadata={}), Star(ra_deg=330.66123248778666, dec_deg=42.255288757286735, metadata={}), Star(ra_deg=330.7162914106881, dec_deg=42.29003046947319, metadata={}), Star(ra_deg=330.7504997672297, dec_deg=42.240428390242094, metadata={}), Star(ra_deg=330.604672301296, dec_deg=42.21859847011704, metadata={}), Star(ra_deg=330.7498637867279, dec_deg=42.28051695514972, metadata={}), Star(ra_deg=330.71430832270755, dec_deg=42.231749645835826, metadata={}), Star(ra_deg=330.6896440316567, dec_deg=42.26675055503721, metadata={}), Star(ra_deg=330.69839011155733, dec_deg=42.233231183945605, metadata={}), Star(ra_deg=330.69230234161955, dec_deg=42.31364995597706, metadata={}), Star(ra_deg=330.65352966410796, dec_deg=42.30382865344687, metadata={}), Star(ra_deg=330.6203455506848, dec_deg=42.27947190278116, metadata={}), Star(ra_deg=330.63271309759074, dec_deg=42.221508687618226, metadata={}), Star(ra_deg=330.7279542229943, dec_deg=42.26749903183207, metadata={}), Star(ra_deg=330.6610720894503, dec_deg=42.28549608518103, metadata={}), Star(ra_deg=330.75031858887354, dec_deg=42.22447686647007, metadata={}), Star(ra_deg=330.6214356288067, dec_deg=42.22166382324433, metadata={}), Star(ra_deg=330.71251399798524, dec_deg=42.24296825730243, metadata={}), Star(ra_deg=330.67146687312385, dec_deg=42.25858000894552, metadata={}), Star(ra_deg=330.68715859371997, dec_deg=42.22927098043517, metadata={}), Star(ra_deg=330.62497084640074, dec_deg=42.30148809133536, metadata={}), Star(ra_deg=330.6379144975921, dec_deg=42.23085018197204, metadata={}), Star(ra_deg=330.6000598354248, dec_deg=42.28111813857364, metadata={}), Star(ra_deg=330.63779624936694, dec_deg=42.2683164896908, metadata={}), Star(ra_deg=330.6956197954143, dec_deg=42.26449713830803, metadata={}), Star(ra_deg=330.6687874155316, dec_deg=42.296568331754905, metadata={}), Star(ra_deg=330.7435608079381, dec_deg=42.30719109560928, metadata={}), Star(ra_deg=330.7437571325056, dec_deg=42.246033893194166, metadata={}), Star(ra_deg=330.750224102727, dec_deg=42.28095636554002, metadata={}), Star(ra_deg=330.6877983773921, dec_deg=42.2455798049503, metadata={}), Star(ra_deg=330.69814907831307, dec_deg=42.21059668481071, metadata={}), Star(ra_deg=330.63296860705526, dec_deg=42.20781412028666, metadata={}), Star(ra_deg=330.63451311645355, dec_deg=42.251433573347576, metadata={}), Star(ra_deg=330.602982083104, dec_deg=42.253868015717075, metadata={}), Star(ra_deg=330.63946346033435, dec_deg=42.28459552979875, metadata={}), Star(ra_deg=330.7270703392626, dec_deg=42.28809688032638, metadata={}), Star(ra_deg=330.6315676630352, dec_deg=42.31652054483156, metadata={}), Star(ra_deg=330.65949552053, dec_deg=42.28831792797967, metadata={}), Star(ra_deg=330.6979625076869, dec_deg=42.31133478915652, metadata={}), Star(ra_deg=330.7183754234475, dec_deg=42.22426159098834, metadata={}), Star(ra_deg=330.68914321047885, dec_deg=42.238512322362915, metadata={}), Star(ra_deg=330.718978380086, dec_deg=42.218431057161126, metadata={}), Star(ra_deg=330.6072142255936, dec_deg=42.29492068220128, metadata={}), Star(ra_deg=330.64958516083357, dec_deg=42.24314816918099, metadata={}), Star(ra_deg=330.6194939327094, dec_deg=42.27653138861859, metadata={}), Star(ra_deg=330.6354925008595, dec_deg=42.26213049391419, metadata={}), Star(ra_deg=330.6546085584999, dec_deg=42.30462078903871, metadata={}), Star(ra_deg=330.6395134548358, dec_deg=42.27428230972779, metadata={}), Star(ra_deg=330.7066908506129, dec_deg=42.25041445019963, metadata={}), Star(ra_deg=330.72020664021653, dec_deg=42.31989124275012, metadata={}), Star(ra_deg=330.6941742605178, dec_deg=42.307609370648855, metadata={}), Star(ra_deg=330.7235707472814, dec_deg=42.256392562999395, metadata={}), Star(ra_deg=330.6951164979502, dec_deg=42.246608653434905, metadata={}), Star(ra_deg=330.69294766722993, dec_deg=42.22849553693228, metadata={}), Star(ra_deg=330.59866319121795, dec_deg=42.2588853079003, metadata={}), Star(ra_deg=330.64790624735826, dec_deg=42.21852811122868, metadata={}), Star(ra_deg=330.63555748468946, dec_deg=42.23690247822715, metadata={}), Star(ra_deg=330.64298653047155, dec_deg=42.29498870969941, metadata={}), Star(ra_deg=330.65091761221714, dec_deg=42.253832762527445, metadata={}), Star(ra_deg=330.6083325601781, dec_deg=42.219466164132726, metadata={}), Star(ra_deg=330.69246722434815, dec_deg=42.30032181091403, metadata={}), Star(ra_deg=330.7047567858923, dec_deg=42.26530788748162, metadata={}), Star(ra_deg=330.73407391006526, dec_deg=42.26570053927472, metadata={}), Star(ra_deg=330.66692426152565, dec_deg=42.28604080283418, metadata={}), Star(ra_deg=330.63227937805783, dec_deg=42.2699770726146, metadata={}), Star(ra_deg=330.66263616367706, dec_deg=42.255171424318846, metadata={}), Star(ra_deg=330.6930368863675, dec_deg=42.28229574770448, metadata={})), quad_stars=(Star(ra_deg=330.6312907472159, dec_deg=42.27082138820124, metadata={}), Star(ra_deg=330.6930368863675, dec_deg=42.28229574770448, metadata={}), Star(ra_deg=330.66123248778666, dec_deg=42.255288757286735, metadata={}), Star(ra_deg=330.66692426152565, dec_deg=42.28604080283418, metadata={})), wcs_fields={'WCSAXES': (2, 'Number of coordinate axes'), 'EQUINOX': (2000.0, 'Equatorial coordinates definition (yr)'), 'LONPOLE': (180.0, 'Native longitude of celestial pole (deg)'), 'LATPOLE': (0.0, 'Native latitude of celestial pole (deg)'), 'CRVAL1': (330.66371747201305, 'RA of reference point'), 'CRVAL2': (42.273373277394306, 'DEC of reference point'), 'CRPIX1': (443.6286307242899, 'X reference pixel'), 'CRPIX2': (338.72823383463896, 'Y reference pixel'), 'CUNIT1': ('deg', 'X pixel scale units'), 'CUNIT2': ('deg', 'Y pixel scale units'), 'CD1_1': (-6.99512038152606e-06, 'Transformation matrix'), 'CD1_2': (0.00014541064154927877, 'Transformation matrix'), 'CD2_1': (0.00014541064154927877, 'Transformation matrix'), 'CD2_2': (6.995120381526028e-06, 'Transformation matrix'), 'CTYPE1': ('RA---TAN-SIP', 'TAN (gnomonic) projection + SIP distortions'), 'CTYPE2': ('DEC--TAN-SIP', 'TAN (gnomonic) projection + SIP distortions')})])

We can see that the solver was successful, and that the result has a logodds much higher than the default minimum of 21, meaning that the calibration was easy and trustable.

The WCS (for one set of pairs) can be obtained with

from astropy.wcs import WCS
wcs1 = WCS(solution.best_match().wcs_fields)

We can plot the image with the sources and the WCS overlayed, and the catalog sources overplotted

fig, ax = plt.subplots(subplot_kw={'projection': wcs1})

imshow_w_sources(redf.mdata, ax=ax)

from iop4lib.db import AstroSource
from photutils.aperture import CircularAperture

for src in AstroSource.get_sources_in_field(wcs1, width=redf.mdata.shape[1], height=redf.mdata.shape[0]):
    ax.plot(*src.coord.to_pixel(wcs1), 'rx')
    ax.annotate(src.name, src.coord.to_pixel(wcs1), xytext=(20,0), textcoords='offset points', color='red', weight='bold')

ax.coords.grid(color='white', linestyle='solid')

ax.coords['ra'].set_axislabel('Right Ascension')
ax.coords['ra'].set_ticklabel_visible(True)
ax.coords['ra'].set_ticklabel_position('lb')

ax.coords['dec'].set_axislabel('Declination')
ax.coords['dec'].set_ticklabel_visible(True)
ax.coords['dec'].set_ticklabel_position('lb')

plt.show()

We can see that the sources are well aligned with the catalog sources. Since we already know the separation between the pairs, we can directly build the other wcs by displacing the first one. However this will be done automatically when we call the higher level methods in the ReducedFit class. In fact, the ReducedFit.astrometric_calibration() method will do all the above and save both WCS in the FITS header of the ReducedFit object:

redf.astrometric_calibration()

2024-09-24 16:03 - pid 4128 - [astrometry.py:179] - DEBUG - <ReducedFit 151 | CAHA-T220/2022-09-18/caf-20220918-23:01:33-sci-agui.fits>: 27 different combinations of parameters to try.
2024-09-24 16:03 - pid 4128 - [astrometry.py:182] - DEBUG - <ReducedFit 151 | CAHA-T220/2022-09-18/caf-20220918-23:01:33-sci-agui.fits>: attempt 1 / 27, ({'keep_n_seg': 300, 'border_margin_px': 5, 'output_logodds_threshold': 14, 'n_rms_seg': 1.5, 'bkg_filter_size': 11, 'bkg_box_size': 16, 'seg_fwhm': 1.0, 'seg_kernel_size': None, 'npixels': 32, 'allsky': False, 'd_eps': 1.0, 'dx_eps': 1.0, 'dy_eps': 1.0, 'dx_min': 29.00547116, 'dx_max': 42.44437116, 'dy_min': -4.8838792, 'dy_max': 5.278269900000001, 'd_min': 30.670903331906715, 'd_max': 40.78002746382967, 'bins': 400, 'hist_range': (0, 400), 'position_hint': PositionHint(ra_deg=330.9182, dec_deg=42.38482, radius_deg=0.4166666666666667)}) ...
2024-09-24 16:03 - pid 4128 - [astrometry.py:275] - DEBUG - <ReducedFit 151 | CAHA-T220/2022-09-18/caf-20220918-23:01:33-sci-agui.fits>: len(pos_seg)=177
2024-09-24 16:03 - pid 4128 - [astrometry.py:282] - DEBUG - <ReducedFit 151 | CAHA-T220/2022-09-18/caf-20220918-23:01:33-sci-agui.fits>: Removing segments within 5 px from border.
2024-09-24 16:03 - pid 4128 - [astrometry.py:289] - DEBUG - <ReducedFit 151 | CAHA-T220/2022-09-18/caf-20220918-23:01:33-sci-agui.fits>: seg pairs -> 91 (51.4%), seg_disp_sign=[-34.36083418  -2.02354451]
2024-09-24 16:03 - pid 4128 - [astrometry.py:291] - DEBUG - <ReducedFit 151 | CAHA-T220/2022-09-18/caf-20220918-23:01:33-sci-agui.fits>: seg pairs best -> 87 (49.2%), seg_disp_sign_best=[-34.84111412  -2.2466433 ]
2024-09-24 16:03 - pid 4128 - [astrometry.py:293] - DEBUG - <ReducedFit 151 | CAHA-T220/2022-09-18/caf-20220918-23:01:33-sci-agui.fits>: seg pairs xy -> 74, disp_sign_xy=[-36.34098767  -0.1852933 ]
2024-09-24 16:03 - pid 4128 - [astrometry.py:295] - DEBUG - <ReducedFit 151 | CAHA-T220/2022-09-18/caf-20220918-23:01:33-sci-agui.fits>: seg pairs xy best -> 74 (41.8%), seg_disp_sign_xy_best=[-36.34098767  -0.1852933 ]
2024-09-24 16:03 - pid 4128 - [astrometry.py:312] - DEBUG - <ReducedFit 151 | CAHA-T220/2022-09-18/caf-20220918-23:01:33-sci-agui.fits>: trying to solve astrometry with Seg Best XY Pairs (n=74) (output_logodds_threshold=14).
2024-09-24 16:03 - pid 4222 - [__init__.py:212] - INFO - loaded 240 index files
2024-09-24 16:03 - pid 4222 - [__init__.py:276] - INFO - solve 1: start
2024-09-24 16:03 - pid 4222 - [__init__.py:276] - INFO - solve 1: slice=[0, 25[ (1 / 2), index="5200/index-5200-14.fits" (1 / 5)
2024-09-24 16:03 - pid 4222 - [__init__.py:276] - INFO - solve 1: logodds=291.605, matches=64, conflicts=0, distractors=5, ra=330.676, dec=42.2691, scale=0.526041, index="5200/index-5200-14.fits"
2024-09-24 16:03 - pid 4128 - [astrometry.py:327] - DEBUG - <ReducedFit 151 | CAHA-T220/2022-09-18/caf-20220918-23:01:33-sci-agui.fits>: Seg Best XY Pairs (n=74) worked.
2024-09-24 16:03 - pid 4128 - [astrometry.py:328] - DEBUG - <ReducedFit 151 | CAHA-T220/2022-09-18/caf-20220918-23:01:33-sci-agui.fits>: bm.index_path=PosixPath('/mnt/astrometry_cache/5200/index-5200-14.fits')
2024-09-24 16:03 - pid 4128 - [astrometry.py:329] - DEBUG - <ReducedFit 151 | CAHA-T220/2022-09-18/caf-20220918-23:01:33-sci-agui.fits>: bm.center_ra_deg=330.6759838494456
2024-09-24 16:03 - pid 4128 - [astrometry.py:330] - DEBUG - <ReducedFit 151 | CAHA-T220/2022-09-18/caf-20220918-23:01:33-sci-agui.fits>: bm.center_dec_deg=42.269073080142824
2024-09-24 16:03 - pid 4128 - [astrometry.py:331] - DEBUG - <ReducedFit 151 | CAHA-T220/2022-09-18/caf-20220918-23:01:33-sci-agui.fits>: bm.scale_arcsec_per_pixel=0.5260410872648312
2024-09-24 16:03 - pid 4128 - [astrometry.py:332] - DEBUG - <ReducedFit 151 | CAHA-T220/2022-09-18/caf-20220918-23:01:33-sci-agui.fits>: bm.logodds=291.6054382324219
2024-09-24 16:03 - pid 4128 - [astrometry.py:192] - DEBUG - <ReducedFit 151 | CAHA-T220/2022-09-18/caf-20220918-23:01:33-sci-agui.fits>: WCS built with attempt 1 / 27 ({'keep_n_seg': 300, 'border_margin_px': 5, 'output_logodds_threshold': 14, 'n_rms_seg': 1.5, 'bkg_filter_size': 11, 'bkg_box_size': 16, 'seg_fwhm': 1.0, 'seg_kernel_size': None, 'npixels': 32, 'allsky': False, 'd_eps': 1.0, 'dx_eps': 1.0, 'dy_eps': 1.0, 'dx_min': 29.00547116, 'dx_max': 42.44437116, 'dy_min': -4.8838792, 'dy_max': 5.278269900000001, 'd_min': 30.670903331906715, 'd_max': 40.78002746382967, 'bins': 400, 'hist_range': (0, 400), 'position_hint': PositionHint(ra_deg=330.9182, dec_deg=42.38482, radius_deg=0.4166666666666667)}).
2024-09-24 16:03 - pid 4128 - [astrometry.py:204] - DEBUG - <ReducedFit 151 | CAHA-T220/2022-09-18/caf-20220918-23:01:33-sci-agui.fits>: building summary images.
2024-09-24 16:03 - pid 4128 - [plotting.py:476] - DEBUG - <ReducedFit 151 | CAHA-T220/2022-09-18/caf-20220918-23:01:33-sci-agui.fits>: plotting astrometry summary image of background substraction results
2024-09-24 16:03 - pid 4128 - [plotting.py:488] - DEBUG - <ReducedFit 151 | CAHA-T220/2022-09-18/caf-20220918-23:01:33-sci-agui.fits>: plotting astrometry summary image of segmentation results
2024-09-24 16:03 - pid 4128 - [plotting.py:582] - DEBUG - <ReducedFit 151 | CAHA-T220/2022-09-18/caf-20220918-23:01:33-sci-agui.fits>: plotting astrometry summary image of astrometry results
2024-09-24 16:03 - pid 4128 - [plotting.py:214] - DEBUG - <ReducedFit 151 | CAHA-T220/2022-09-18/caf-20220918-23:01:33-sci-agui.fits>: found 5 catalog sources in field: [AstroSource.objects.get(name='2200+420'), AstroSource.objects.get(name='BL Lacertae C'), AstroSource.objects.get(name='BL Lacertae H'), AstroSource.objects.get(name='BL Lacertae test 2'), AstroSource.objects.get(name='BL Lacertae test 4')]
2024-09-24 16:03 - pid 4128 - [connectionpool.py:1051] - DEBUG - Starting new HTTPS connection (1): simbad.cds.unistra.fr:443
2024-09-24 16:03 - pid 4128 - [connectionpool.py:546] - DEBUG - https://simbad.cds.unistra.fr:443 "POST /simbad/sim-script HTTP/11" 200 None
2024-09-24 16:03 - pid 4128 - [instrument.py:425] - DEBUG - <ReducedFit 151 | CAHA-T220/2022-09-18/caf-20220918-23:01:33-sci-agui.fits>: saving WCSs to FITS header.

redf.wcs1

WCS Keywords

Number of WCS axes: 2
CTYPE : 'RA---TAN-SIP'  'DEC--TAN-SIP'  
CRVAL : 330.71400656577  42.244431367274  
CRPIX : 236.6813171671  588.25901065941  
PC1_1 PC1_2  : -4.7575342185789e-06  0.00014604505454994  
PC2_1 PC2_2  : 0.00014604505454994  4.757534218579e-06  
CDELT : 1.0  1.0  
NAXIS : 800  800

redf.wcs2

WCS Keywords

Number of WCS axes: 2
CTYPE : 'RA---TAN-SIP'  'DEC--TAN-SIP'  
CRVAL : 330.71400656577  42.244431367274  
CRPIX : 200.34032949831  588.07371735926  
PC1_1 PC1_2  : -4.7575342185789e-06  0.00014604505454994  
PC2_1 PC2_2  : 0.00014604505454994  4.757534218579e-06  
CDELT : 1.0  1.0  
NAXIS : 800  800

fig, ax = plt.subplots(subplot_kw={'projection': wcs1})

imshow_w_sources(redf.mdata, ax=ax)

from iop4lib.db import AstroSource
from photutils.aperture import CircularAperture

for src in AstroSource.get_sources_in_field(wcs1, width=redf.mdata.shape[1], height=redf.mdata.shape[0]):
    CircularAperture([*src.coord.to_pixel(redf.wcs1)], r=20).plot(ax=ax, color='r')
    ax.annotate(src.name, src.coord.to_pixel(redf.wcs1), xytext=(20,0), textcoords='offset points', color='red', weight='bold', verticalalignment='center')
    ax.plot(*src.coord.to_pixel(redf.wcs2), 'rx')

ax.coords.grid(color='white', linestyle='solid')

ax.coords['ra'].set_axislabel('Right Ascension')
ax.coords['ra'].set_ticklabel_visible(True)
ax.coords['ra'].set_ticklabel_position('bl')

ax.coords['dec'].set_axislabel('Declination')
ax.coords['dec'].set_ticklabel_visible(True)
ax.coords['dec'].set_ticklabel_position('lb')
#ax.coords['dec'].set_ticklabel(rotation='vertical')

plt.show()

If the calibration was successful, summary images containing more or less the same information as the ones we plotted will be created and should be displayed in the admin page for the reduced fit:

/iop4/admin/iop4api/reducedfit/details/151