使用 astropy 和 fit_wcs_from_points 拟合天顶等面积投影

如何解决使用 astropy 和 fit_wcs_from_points 拟合天顶等面积投影

我正在尝试使用 astropy.wcs.utils.fit_wcs_from_points 来拟合使用天顶等面积投影（WCS 代码 ZEA）投影的点。这种投影在全天空相机中很受欢迎。

我首先用已知的 WCS 投影一组天体坐标，看看我是否可以恢复它。输入值通过以下方式获得：

x,y = w.world_to_pixel(lon * u.deg,lat * u.deg)
world_coords = SkyCoord(lon * u.deg,lat * u.deg)

投影是：

from astropy import wcs
w = wcs.WCS(naxis=2)
scale = 0.095
w.wcs.crpix = [1290,1950]
w.wcs.cdelt = [scale,scale]
w.wcs.crval = [0,90]
w.wcs.ctype = ["ALON-ZEA","ALAT-ZEA"]

在我所有的测试中，当我执行拟合时，我得到以下异常：

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-28-cfb0b01c0d18> in <module>
----> 1 astropy.wcs.utils.fit_wcs_from_points([x,y],world_coords,2                                       proj_point=SkyCoord(0 * u.deg,90 * u.deg),3                                       projection='ZEA')

/usr/lib64/python3.9/site-packages/astropy/wcs/utils.py in fit_wcs_from_points(xy,proj_point,projection,sip_degree)
   1076     # and cd terms are way off.
   1077     p0 = np.concatenate([wcs.wcs.cd.flatten(),wcs.wcs.crpix.flatten()])
-> 1078     fit = least_squares(_linear_wcs_fit,p0,1079                         args=(lon,lat,xp,yp,wcs))
   1080     wcs.wcs.crpix = np.array(fit.x[4:6])

/usr/lib64/python3.9/site-packages/scipy/optimize/_lsq/least_squares.py in least_squares(fun,x0,jac,bounds,method,ftol,xtol,gtol,x_scale,loss,f_scale,diff_step,tr_solver,tr_options,jac_sparsity,max_nfev,verbose,args,kwargs)
    812 
    813     if not np.all(np.isfinite(f0)):
--> 814         raise ValueError("Residuals are not finite in the initial point.")
    815 
    816     n = x0.size

ValueError: Residuals are not finite in the initial point.

我尝试过的东西，没有改变：

通过 projection 而不是 ZEA 中的实际测试 WCS 转换
在我的测试 WCS 中将 CRPIX 设置为 0，因此所有点都在 (0,0) 附近
删除proj_point
尝试使用 astropy 4.2.1（最新发布）

此示例代码重现了该问题：

import numpy as np
import astropy.wcs.utils
from astropy.coordinates import SkyCoord
import astropy.units as u

x0 = np.array([702.4,1480.4,1223.5,897,1916.6])
y0 = np.array([1925.8,2269.3,2679.1,1632.7,1586.3])
zea_lon = [268.8,145.6,181.6,305.4,56.5]
zea_lat = [31.8,54.0,15.2,40.6,16.1]
world_coords0 = SkyCoord(zea_lon * u.deg,zea_lat * u.deg)
astropy.wcs.utils.fit_wcs_from_points([x0,y0],world_coords0,proj_point=SkyCoord(0 * u.deg,projection='ZEA')