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()
2026-06-23 15:23 - pid 3865 - [instrument.py:540] - DEBUG - <ReducedFit 118 | OSN-T090/2022-09-18/BLLac-0001R.fit>: building file
2026-06-23 15:23 - pid 3865 - [instrument.py:446] - DEBUG - <ReducedFit 118 | OSN-T090/2022-09-18/BLLac-0001R.fit>: applying masters
2026-06-23 15:23 - pid 3865 - [instrument.py:455] - WARNING - <ReducedFit 118 | OSN-T090/2022-09-18/BLLac-0001R.fit>: no masterdark found, assuming dark current = 0, is this a CCD camera and it's cold?
2026-06-23 15:23 - pid 3865 - [instrument.py:546] - DEBUG - <ReducedFit 118 | OSN-T090/2022-09-18/BLLac-0001R.fit>: performing astrometric calibration
2026-06-23 15:23 - pid 3865 - [astrometry.py:214] - DEBUG - <ReducedFit 118 | OSN-T090/2022-09-18/BLLac-0001R.fit>: 27 different combinations of parameters to try.
2026-06-23 15:23 - pid 3865 - [astrometry.py:217] - DEBUG - <ReducedFit 118 | OSN-T090/2022-09-18/BLLac-0001R.fit>: attempt 1 / 27, ({'position_hint': PositionHint(ra_deg=330.7291666666667, dec_deg=42.26611111111111, radius_deg=0.32999999999999996), '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}) ...
2026-06-23 15:23 - pid 3865 - [astrometry.py:299] - DEBUG - <ReducedFit 118 | OSN-T090/2022-09-18/BLLac-0001R.fit>: len(pos_seg)=187
2026-06-23 15:23 - pid 3865 - [astrometry.py:302] - DEBUG - <ReducedFit 118 | OSN-T090/2022-09-18/BLLac-0001R.fit>: Keeping only 150 brightest segments.
2026-06-23 15:23 - pid 3865 - [astrometry.py:306] - DEBUG - <ReducedFit 118 | OSN-T090/2022-09-18/BLLac-0001R.fit>: Removing segments within 5 px from border.
2026-06-23 15:23 - pid 3865 - [astrometry.py:349] - DEBUG - <ReducedFit 118 | OSN-T090/2022-09-18/BLLac-0001R.fit>: trying to solve astrometry with Seg Pos (n=147) (output_logodds_threshold=14).
2026-06-23 15:23 - pid 3894 - [__init__.py:215] - INFO - loaded 240 index files
2026-06-23 15:23 - pid 3894 - [__init__.py:297] - INFO - solve 1: start
2026-06-23 15:23 - pid 3894 - [__init__.py:297] - INFO - solve 1: slice=[0, 25) (1 / 5), index="5200/index-5200-14.fits" (1 / 5)
2026-06-23 15:23 - pid 3894 - [__init__.py:297] - 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"
2026-06-23 15:23 - pid 3865 - [astrometry.py:364] - DEBUG - <ReducedFit 118 | OSN-T090/2022-09-18/BLLac-0001R.fit>: Seg Pos (n=147) worked.
2026-06-23 15:23 - pid 3865 - [astrometry.py:365] - DEBUG - <ReducedFit 118 | OSN-T090/2022-09-18/BLLac-0001R.fit>: bm.index_path=PosixPath('/mnt/astrometry_cache/5200/index-5200-14.fits')
2026-06-23 15:23 - pid 3865 - [astrometry.py:366] - DEBUG - <ReducedFit 118 | OSN-T090/2022-09-18/BLLac-0001R.fit>: bm.center_ra_deg=330.6954314398751
2026-06-23 15:23 - pid 3865 - [astrometry.py:367] - DEBUG - <ReducedFit 118 | OSN-T090/2022-09-18/BLLac-0001R.fit>: bm.center_dec_deg=42.26972920042598
2026-06-23 15:23 - pid 3865 - [astrometry.py:368] - DEBUG - <ReducedFit 118 | OSN-T090/2022-09-18/BLLac-0001R.fit>: bm.scale_arcsec_per_pixel=0.7736712092056637
2026-06-23 15:23 - pid 3865 - [astrometry.py:369] - DEBUG - <ReducedFit 118 | OSN-T090/2022-09-18/BLLac-0001R.fit>: bm.logodds=380.9163513183594
2026-06-23 15:23 - pid 3865 - [astrometry.py:222] - DEBUG - <ReducedFit 118 | OSN-T090/2022-09-18/BLLac-0001R.fit>: WCS built with attempt 1 / 27 ({'position_hint': PositionHint(ra_deg=330.7291666666667, dec_deg=42.26611111111111, radius_deg=0.32999999999999996), '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}).
2026-06-23 15:23 - pid 3865 - [astrometry.py:234] - DEBUG - <ReducedFit 118 | OSN-T090/2022-09-18/BLLac-0001R.fit>: building summary images.
2026-06-23 15:23 - pid 3865 - [plotting.py:475] - DEBUG - <ReducedFit 118 | OSN-T090/2022-09-18/BLLac-0001R.fit>: plotting astrometry summary image of background substraction results
2026-06-23 15:23 - pid 3865 - [plotting.py:487] - DEBUG - <ReducedFit 118 | OSN-T090/2022-09-18/BLLac-0001R.fit>: plotting astrometry summary image of segmentation results
2026-06-23 15:23 - pid 3865 - [plotting.py:587] - DEBUG - <ReducedFit 118 | OSN-T090/2022-09-18/BLLac-0001R.fit>: plotting astrometry summary image of astrometry results
2026-06-23 15:23 - pid 3865 - [plotting.py:213] - DEBUG - <ReducedFit 118 | OSN-T090/2022-09-18/BLLac-0001R.fit>: found 31 catalog sources in field: <AstroSourceQuerySet [AstroSource.objects.get(name='2200+420'), AstroSource.objects.get(name='PanSTARRS 330.70127 42.31286'), AstroSource.objects.get(name='PanSTARRS 330.65424 42.22716'), AstroSource.objects.get(name='PanSTARRS 330.69304 42.28230'), AstroSource.objects.get(name='PanSTARRS 330.60371 42.27227'), AstroSource.objects.get(name='PanSTARRS 330.78031 42.29269'), AstroSource.objects.get(name='PanSTARRS 330.74253 42.30545'), AstroSource.objects.get(name='PanSTARRS 330.68965 42.26675'), AstroSource.objects.get(name='PanSTARRS 330.77540 42.22265'), AstroSource.objects.get(name='PanSTARRS 330.61553 42.24694'), AstroSource.objects.get(name='PanSTARRS 330.71458 42.25859'), AstroSource.objects.get(name='PanSTARRS 330.65030 42.28163'), AstroSource.objects.get(name='PanSTARRS 330.67245 42.31685'), AstroSource.objects.get(name='PanSTARRS 330.67544 42.21823'), AstroSource.objects.get(name='PanSTARRS 330.63461 42.19721'), AstroSource.objects.get(name='PanSTARRS 330.67445 42.21723'), AstroSource.objects.get(name='PanSTARRS 330.60190 42.33041'), AstroSource.objects.get(name='PanSTARRS 330.63916 42.24036'), AstroSource.objects.get(name='PanSTARRS 330.58896 42.22607'), AstroSource.objects.get(name='PanSTARRS 330.62895 42.20273'), '...(remaining elements truncated)...']>
2026-06-23 15:23 - pid 3865 - [connectionpool.py:1062] - DEBUG - Starting new HTTPS connection (1): simbad.cds.unistra.fr:443
2026-06-23 15:23 - pid 3865 - [connectionpool.py:544] - DEBUG - https://simbad.cds.unistra.fr:443 "GET /simbad/sim-tap/capabilities HTTP/1.1" 200 None
2026-06-23 15:23 - pid 3865 - [connectionpool.py:544] - DEBUG - https://simbad.cds.unistra.fr:443 "POST /simbad/sim-tap/sync HTTP/1.1" 200 None
2026-06-23 15:23 - pid 3865 - [connectionpool.py:544] - DEBUG - https://simbad.cds.unistra.fr:443 "POST /simbad/sim-tap/sync HTTP/1.1" 200 None
2026-06-23 15:23 - pid 3865 - [connectionpool.py:544] - DEBUG - https://simbad.cds.unistra.fr:443 "POST /simbad/sim-tap/sync HTTP/1.1" 200 None
2026-06-23 15:23 - pid 3865 - [logger.py:235] - WARNING - AstropyDeprecationWarning: "epoch" was deprecated in version 0.4.8 and will be removed in a future version.
2026-06-23 15:23 - pid 3865 - [connectionpool.py:544] - DEBUG - https://simbad.cds.unistra.fr:443 "POST /simbad/sim-tap/sync HTTP/1.1" 200 None
2026-06-23 15:23 - pid 3865 - [instrument.py:485] - DEBUG - <ReducedFit 118 | OSN-T090/2022-09-18/BLLac-0001R.fit>: saving WCSs to FITS header.
2026-06-23 15:23 - pid 3865 - [instrument.py:557] - DEBUG - <ReducedFit 118 | OSN-T090/2022-09-18/BLLac-0001R.fit>: astrometric calibration was successful.
2026-06-23 15:23 - pid 3865 - [instrument.py:560] - DEBUG - <ReducedFit 118 | OSN-T090/2022-09-18/BLLac-0001R.fit>: searching for sources in field...
2026-06-23 15:23 - pid 3865 - [instrument.py:563] - DEBUG - <ReducedFit 118 | OSN-T090/2022-09-18/BLLac-0001R.fit>: found 31 sources in field.
WARNING: AstropyDeprecationWarning: "epoch" was deprecated in version 0.4.8 and will be removed in a future version. [iop4lib.utils]
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
- telescope: OSN-T090
- night: 2022-09-18
- filename: BLLac-0001R.fit
- instrument: AndorT90
- imgtype: LIGHT
- imgbinning: 2x2
- 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")
---------------------------------------------------------------------------
DoesNotExist Traceback (most recent call last)
Cell In[6], line 2
1 from iop4lib.db import RawFit
----> 2 rf = RawFit.by_fileloc("CAHA-T220/2022-09-18/caf-20220918-23:01:33-sci-agui.fits")
File ~/actions-runner/_work/iop4/iop4/iop4lib/db/rawfit.py:370, in RawFit.by_fileloc(cls, fileloc)
367 @classmethod
368 def by_fileloc(cls, fileloc):
369 telescope, night, filename = cls.fileloc_to_tel_night_filename(fileloc)
--> 370 return cls.objects.filter(epoch__telescope=telescope, epoch__night=night, filename=filename).get()
File ~/actions-runner/_work/_tool/Python/3.11.15/x64/lib/python3.11/site-packages/django/db/models/query.py:639, in QuerySet.get(self, *args, **kwargs)
637 return clone._result_cache[0]
638 if not num:
--> 639 raise self.model.DoesNotExist(
640 "%s matching query does not exist." % self.model._meta.object_name
641 )
642 raise self.model.MultipleObjectsReturned(
643 "get() returned more than one %s -- it returned %s!"
644 % (
(...) 647 )
648 )
DoesNotExist: RawFit matching query does not exist.
We first create the reduced image and apply the calibration frames
redf = ReducedFit.create(rawfit=rf)
redf.apply_masters()
The .has_pairs property tells whether an image is expected to have pairs:
redf.has_pairs
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
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
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
redf.header_hintobject.coord.separation(redf.header_hintcoord)
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
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
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
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)
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)
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
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(wcs=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()
redf.wcs1
redf.wcs2
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(wcs=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: