I'm trying to use astropy.wcs.utils.fit_wcs_from_points to fit points projected with zenithal equal area projection (WCS code ZEA). This projection is popular in all-sky cameras.
I have started by projecting a set of celestial coordinates with a known WCS, to see if I can recover it. Input values are obtained by:
x, y = w.world_to_pixel(lon * u.deg, lat * u.deg)
world_coords = SkyCoord(lon * u.deg, lat * u.deg)
and the projection is:
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"]
In all my tests I get the following exception when I perform the fitting:
---------------------------------------------------------------------------
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, world_coords, 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.
Things I have tried, without changes:
- pass the actual test WCS transformation in
projection, instead ofZEA - set
CRPIXto 0 in my test WCS, so all the points are around (0 ,0) - remove
proj_point - try with astropy 4.2.1 (latest released)
This sample code reproduces the problem:
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, 90 * u.deg),
projection='ZEA')