Skip to content

Automated distortion fits using WCS file provided by astrometry.net #535

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
dvida wants to merge 2 commits into from
Closed

Automated distortion fits using WCS file provided by astrometry.net #535

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
890d9df
Began work implementing automated distortion fits using WCS file prov…
dvida Mar 12, 2025
08bcdba
fix: update polynomial distortion order to 3 in novaAstrometryNetSolve
dvida Mar 12, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 7 additions & 1 deletion RMS/Astrometry/AstrometryNet.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -242,7 +242,13 @@ def astrometryNetSolve(ff_file_path=None, img=None, mask=None, x_data=None, y_da
fov_h = y_max*solution.best_match().scale_arcsec_per_pixel/3600


star_data = [x_data, y_data]
# Compute the RA and Dec of x_data, y_data
x_data = np.array(x_data)
y_data = np.array(y_data)

ra_data, dec_data = wcs_obj.all_pix2world(x_data, y_data, 1)

star_data = [x_data, y_data, ra_data, dec_data]

return ra_mid, dec_mid, rot_eq_standard, scale, fov_w, fov_h, star_data

Expand Down
45 changes: 44 additions & 1 deletion RMS/Astrometry/AstrometryNetNova.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -368,7 +368,11 @@ def _printWebLink(stat, first_status=None):
# Add keyword arguments
kwargs = {}
kwargs['publicly_visible'] = 'y'

# Force the center to be at the center of the image
kwargs['crpix_center'] = True

# Set the polynomial distortion order to 3
kwargs['tweak_order'] = 3

# Add the scale to keyword arguments, if given
Expand Down Expand Up @@ -544,7 +548,46 @@ def _printWebLink(stat, first_status=None):
fov_w = result['width_arcsec']/3600
fov_h = result['height_arcsec']/3600

return ra_mid, dec_mid, rot_eq_standard, scale, fov_w, fov_h, None

### Construct an array of star data ###

# If the star data is available, use it and map to RA and dec
if (x_data is not None) and (y_data is not None):

# Convert the star data to RA and dec
ra_data, dec_data = wcs_obj.all_pix2world(x_data, y_data, 1)

star_data = [x_data, y_data, ra_data, dec_data]

elif img is not None:

# If the star x and y coordinates are not available, sample image coordinates are regular intervals
# to get the RA and dec coordinates

# Image dimensions
h, w = img.shape

# Start 10 px from the edge, sample 5 samples in each image dimension for a total of 25 samples
x_data = np.linspace(10, w-10, 5)
y_data = np.linspace(10, h-10, 5)

x_data, y_data = np.meshgrid(x_data, y_data)
x_data = x_data.flatten()
y_data = y_data.flatten()


# Convert the star data to RA and dec
ra_data, dec_data = wcs_obj.all_pix2world(x_data, y_data, 1)

star_data = [x_data, y_data, ra_data, dec_data]

else:
star_data = None

###


return ra_mid, dec_mid, rot_eq_standard, scale, fov_w, fov_h, star_data


if __name__ == '__main__':
Expand Down
14 changes: 7 additions & 7 deletions RMS/Formats/Platepar.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -411,7 +411,7 @@ def _calcImageResidualsAstro(params, platepar, jd, catalog_stars, img_stars):
if not fixed_scale:
F_scale = params[3]

img_x, img_y, _ = img_stars.T
img_x, img_y = img_stars[:, :2].T

pp_copy = copy.deepcopy(platepar)

Expand Down Expand Up @@ -441,7 +441,7 @@ def _calcSkyResidualsAstro(params, platepar, jd, catalog_stars, img_stars):
if not fixed_scale:
F_scale = params[3]

img_x, img_y, _ = img_stars.T
img_x, img_y = img_stars[:, :2].T

pp_copy = copy.deepcopy(platepar)

Expand All @@ -453,7 +453,7 @@ def _calcSkyResidualsAstro(params, platepar, jd, catalog_stars, img_stars):
if not fixed_scale:
pp_copy.F_scale = abs(F_scale)

img_x, img_y, _ = img_stars.T
img_x, img_y = img_stars[:, :2].T

# Get image coordinates of catalog stars
ra_array, dec_array = getPairedStarsSkyPositions(img_x, img_y, jd, pp_copy)
Expand Down Expand Up @@ -546,7 +546,7 @@ def _calcImageResidualsDistortion(params, platepar, jd, catalog_stars, img_stars
pp_copy.x_poly_rev = np.zeros(platepar.poly_length)
pp_copy.y_poly_rev = params

img_x, img_y, _ = img_stars.T
img_x, img_y = img_stars[:, :2].T

# Get image coordinates of catalog stars
catalog_x, catalog_y, catalog_mag = getCatalogStarsImagePositions(catalog_stars, jd, pp_copy)
Expand Down Expand Up @@ -587,7 +587,7 @@ def _calcSkyResidualsDistortion(params, platepar, jd, catalog_stars, img_stars,
else:
pp_copy.y_poly_fwd = params

img_x, img_y, _ = img_stars.T
img_x, img_y = img_stars[:, :2].T

# Get image coordinates of catalog stars
ra_array, dec_array = getPairedStarsSkyPositions(img_x, img_y, jd, pp_copy)
Expand Down Expand Up @@ -638,7 +638,7 @@ def _calcImageResidualsAstroAndDistortionRadial(params, platepar, jd, catalog_st
# Assign distortion parameters
pp_copy.x_poly_rev = params[4:]

img_x, img_y, _ = img_stars.T
img_x, img_y = img_stars[:, :2].T

# Get image coordinates of catalog stars
catalog_x, catalog_y, catalog_mag = getCatalogStarsImagePositions(catalog_stars, jd, pp_copy)
Expand Down Expand Up @@ -672,7 +672,7 @@ def _calcSkyResidualsAstroAndDistortionRadial(params, platepar, jd, catalog_star
# Assign distortion parameters
pp_copy.x_poly_fwd = params[4:]

img_x, img_y, _ = img_stars.T
img_x, img_y = img_stars[:, :2].T

# Get image coordinates of catalog stars
ra_array, dec_array = getPairedStarsSkyPositions(img_x, img_y, jd, pp_copy)
Expand Down
101 changes: 80 additions & 21 deletions Utils/SkyFit2.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4356,36 +4356,95 @@ def getInitialParamsAstrometryNet(self, upload_image=True):
print(' FOV = {:.2f} x {:.2f} deg'.format(fov_w, fov_h))


# If a list of detected stars is provided by the astrometry.net, use it to run FFT alignment
# If a list of detected stars is provided by the astrometry.net:
# - use it to run FFT alignment
# - fit a platepar with a radial3-odd distortion model
if star_data is not None:

# Construct an array with star coordinates (x, y, ra, dec per row)
x_data, y_data, ra_data, dec_data = star_data


### TEST

# Print the star data
print()
print("Running FFT alignment...")
print("Star data for automated fitting:")
print(" X Y RA Dec")
for x, y, ra, dec in zip(x_data, y_data, ra_data, dec_data):
print("{:8.2f} {:8.2f} {:8.2f} {:8.2f}".format(x, y, ra, dec))

# # Plot the x, y coordinates
# plt.figure()
# plt.scatter(x_data, y_data, c='r', s=10)
# plt.show()

# # Plot the ra, dec coordinates on a polar plot
# plt.figure()
# plt.polar(np.radians(ra_data), 90-dec_data, 'ro')
# plt.show()

# Construct an array with star coordinates (x, y per row)
calstars_coords = np.array(star_data).T

###

calstars_coords = np.array([x_data, y_data]).T
catalog_coords = np.array([ra_data, dec_data, np.zeros_like(ra_data)]).T

# Get the time of the image
calstars_time = self.img_handle.currentTime()

self.platepar = alignPlatepar(
self.config, self.platepar,
calstars_time, calstars_coords,
scale_update=True, show_plot=False
)

self.platepar.updateRefRADec(skip_rot_update=True)

# Make sure there are at least 4 stars for FOV alignment
if len(calstars_coords) > 4:

print()
print("Running FFT alignment...")

self.platepar = alignPlatepar(
self.config, self.platepar,
calstars_time, calstars_coords,
scale_update=True, show_plot=False
)

self.platepar.updateRefRADec(skip_rot_update=True)

print()
print('FFT aligned:')
print('------------------------')
print(' RA = {:.2f} deg'.format(self.platepar.RA_d))
print(' Dec = {:.2f} deg'.format(self.platepar.dec_d))
print(' Azim = {:.2f} deg'.format(self.platepar.az_centre))
print(' Alt = {:.2f} deg'.format(self.platepar.alt_centre))
print(' Rot horiz = {:.2f} deg'.format(self.platepar.rotation_from_horiz))
print(' Pos angle = {:.2f} deg'.format(self.platepar.pos_angle_ref))
print(' Scale = {:.3f} arcmin/px'.format(60/self.platepar.F_scale))


# Make sure there are at least 6 stars for distortion fitting
if len(calstars_coords) > 6:

### Fit a platepar with a radial3-odd distortion model ###
print()
print("Fitting a platepar with a radial3-odd distortion model...")

# Turn on asymmetry correction
self.platepar.asymmetry_corr = True

# Switch the platepar to the radial3-odd distortion model
self.platepar.setDistortionType('radial3-odd', reset_params=True)

# Do a platepar fit on the provided star data
self.platepar.fitAstrometry(
date2JD(*calstars_time),
calstars_coords, catalog_coords,
first_platepar_fit=True, fit_only_pointing=False, fixed_scale=False
)

# Print platepar parameters
print()
print('Fitted platepar:')
print(self.platepar)

print()
print('FFT aligned:')
print('------------------------')
print(' RA = {:.2f} deg'.format(self.platepar.RA_d))
print(' Dec = {:.2f} deg'.format(self.platepar.dec_d))
print(' Azim = {:.2f} deg'.format(self.platepar.az_centre))
print(' Alt = {:.2f} deg'.format(self.platepar.alt_centre))
print(' Rot horiz = {:.2f} deg'.format(self.platepar.rotation_from_horiz))
print(' Pos angle = {:.2f} deg'.format(self.platepar.pos_angle_ref))
print(' Scale = {:.3f} arcmin/px'.format(60/self.platepar.F_scale))


def getFOVcentre(self):
Expand Down