Fit_wcs_from_points() function raises ‘ValueError: Residuals are not finite in the initial point.'

Hi all, I’m a miserable undergrad working on my observation course data. I’m trying to use astropy’s fit_wcs_from_points() to produce a wcs for my frame, but received the error ‘ValueError: Residuals are not finite in the initial point.’ I tracked down the error in pdb and source code and found that wcs.wcs.p2s(xy, 0) would produce results with lots of nan. I’ve included a minimal working example and a detailed process of how I tracked down the error in this question.py file.
I have dug up stackoverflow and found only one similar question, but it was left unanswered. The link to this question is also included in question.py. I’ve also tried using different reference points or increasing the number of points used, but it didn’t work either.

Minimal working example:

from astropy.wcs.utils import fit_wcs_from_points
from astropy.coordinates import SkyCoord
c = SkyCoord(["11:46:53.867 +72:21:19.64", "11:46:16.550 +72:15:53.89", "11:43:32.195 +72:04:56.81"],frame = 'icrs', unit=(u.hourangle,u.deg))
xy = (array([1050.93327996, 1682.54476132, 3802.23270946]), array([ 809.96831739, 1345.84119511, 2066.74814281]))
#Known physical coordinates of three targets, manually matched from observation frame and Vizier map
w = fit_wcs_from_points(xy, c, projection = 'SIN')
#Where the value has arised. ValueError: Residuals are not finite in the initial point.

I checked fit_wcs_from_points() with pdb and found that the error arised from fit = least_squares(_linear_wcs_fit,p0,args=(lon, lat, xp, yp, wcs),bounds=[[-np.inf, -np.inf, -np.inf, -np.inf, xpmin + 1, ypmin + 1],[np.inf, np.inf, np.inf, np.inf, xpmax + 1, ypmax + 1]])

In this least_squares() function, the error further rise from _linear_wcs_fit(p0, lon, lat, xp, yp, wcs), which returned array([nan, nan, nan, nan, nan, nan])

Further down, error is produced at lon2, lat2 = wcs.wcs_pix2world(xp, yp, 0), where it returns lat2=[nan, nan, nan] (lon2 should probably be the same though i didin’t check)

Return value of wcs_pix2world() is self._array_converter(lambda xy, o: self.wcs.p2s(xy, o)[‘world’],‘output’, *args, **kwargs), so what really went wrong was self.wcs.p2s(xy, o).

Deepest level minimal Working Example:
Executing w = fit_wcs_from_points in ipython would produce a wcs as the program stops with the error,with the following information
WCS Keywords
#Number of WCS axes: 2
#CTYPE : ‘RA—SIN’ ‘DEC–SIN’
#CRVAL : 176.3074171455784 72.21939932332155
#CRPIX : 2426.5829947141146 1438.3582300998578
#CD1_1 CD1_2 : 1.0 0.0
#CD2_1 CD2_2 : 0.0 1.0
#NAXIS : 3803 2067
with this wcs:

wcs.wcs.p2s(anchors, 0)

would produce {‘imgcrd’: array([[nan, nan, nan],[nan, nan, nan]]), ‘phi’: array([nan, nan]), ‘theta’: array([nan, nan]), ‘world’: array([[nan, nan, 0.],[nan, nan, 0.]]), ‘stat’: array([3, 3], dtype=int32)}

I could not find the definition of p2s() anywhere online, so have to stop here. Though I did find the handbook at Wcsprm — Astropy v6.0.0 But I’m not sure if this is the p2s() used, as I did not see the keyword Wcsprm in the source code

Other things I’ve tried:

  • Using different reference points other than these three
  • Increasing number of reference points up to 22. These data might be found in a second anchor.txt file

A similar question can be found on stackoverflow, though no one answered it python - Using astropy fit_wcs_from_points to give FITS file a new WCS - Stack Overflow

Any help is appreciated. Happy New Year!