Mismatch in how astrometric offsets are calculated and how they are applied

[AkashA98]

In general, if we want the offsets to be added additively, i.e RA_NEW = RA + RAOFF and DEC_NEW = DEC + DECOFF they must be calculated in a global frame (a frame that does not depend on the sky position).

But the way astrometric offsets are estimated in this pipeline is by using spherical_offsets_to (see .coorections.crossmatch_qtables --- here) which estimates lat and lon in a local frame. But they are applied in a purely additive fashion (here), which yields wrong results. Offset corrections need to be transformed back to a global frame if we want them to be purely additive.

which is to say (the way spherical_offsets_to in astropy works)

If T is the transformation matrix that transforms to a local frame (where local frame is defined as a frame in which the given coordinates get local_ra = 0 and local_dec = 0), then the offsets calculated currently are the lat and lon of the target point (Q) in the local frame, i.e dra, ddec = lon(TQ), lat(TQ)

So to apply them, we need to build a pointing vector in this local frame r = spherical(theta=ddec, phi=dra) and then get back Q as matrix_multplication(inverse(T), r)

This is all internally done by spherical_offsets_by (here, but using this astropy function as is means that we can not propagate errors in a straight-forward manner (we need to rely on numerical differentiation of this function to apply the fitting errors).

This can be resolved in two ways:

1. Use analytical expressions involving matrices as described above (might add computational overhead, especially for catalogs where the number of sources can be large)
2. Do the offset calculation in a global frame, not in local frame.

[AkashA98]

Currently, I am inclined to go with 2 and re-write the function that does the offset estimation.

[AkashA98]

75 should fix this.

[dlakaplan]

I think applying the corrections including the cos(Dec) factor should be fine. No rotation/transformation matrix is needed. It's true that this is an approximation when the FOV is very large (e.g., which Dec to use) but I am not convinced that this PR is correct.

[AkashA98]

Ok, @dlakaplan, I think I found out where I was going wrong. Yes, a pure difference (just the difference in RA's) is not correct and we need to evaluate that in the tangent plane, but I figured out the actual correction. If I proceed with all the matrix conversion (which essentially does the same tangent plane calculation), we get the RA correction as

$tan(\Delta RA) = \frac{cos(\phi) * sin(\Delta \lambda)}{cos^2(\phi) * cos(\Delta \lambda) + sin^2(\phi)}$

where $\lambda$ and $\phi$ are longitude and latitude, $\Delta RA$ is the RA offset evaluated in the tangent plane and $\Delta \lambda$ is the actual RA correction

which in the low angle limit $\Delta \lambda << 1$

this reduces to

$\Delta \lambda = \Delta RA /cos(\phi)$ (which is what is done in the code).

But since I did not realize this and started working with large offsets (in degrees), I was getting those discrepancies that we discussed.

I did a comparison test of these two solutions and here is what it looks like (it essentially comes down to what you said -- for arcsec offsets that this code is supposed to work with, it doesn't make any difference, but if someone unknowingly uses it to correct fields with ~ degree corrections, then they would get biased (on arcsec level) results)

Here is for a 1 degree wide field:

1. dec = 15
2. dec = 45
3. dec = 89.5

For a 5 degree wide field

1. dec = 15
2. dec=45
3. dec=89.5

So, yes, as you said, for arcsec corrections, the whole matrix story I wrote above (unknowingly) is the same as what the code approximates (with the cos(dec) factor), but since I was testing with large angle offsets, I mistakenly thought what the code did was inconsistent.

I can replace the approximate expression with the fully correct expression, with little effect for arcsec corrections or just put a warning statement with the approximate correction saying that the code is not intended to be used for large corrections.

[AkashA98]

Also, @dlakaplan @ddobie , sorry for the confusion.

[ddobie]

I think either option is probably fine - if you've already written the code to do the full corrections then we might as well use it, but if it's not written I think a warning is sufficient!

[AkashA98]

And here is the code for reference:

`

dec = 45
s1 = SkyCoord(ra=0*u.degree, dec=dec*u.degree)
s2 = SkyCoord(ra=np.arange(0.1, 5 * 3600.1, 1)*u.arcsec, dec=dec*u.degree)
off_ra, off_dec = s1.spherical_offsets_to(s2)
adj_ra = off_ra/np.cos(dec*u.degree)
err1 = (-adj_ra +s2.ra).to(u.arcsec)

# Write the new method
lam = np.tan(off_ra)
epsilon = 1.e-1 * u.arcsec
err = np.abs(lam**2 * np.cos(2*s1.dec) / 4 / np.cos(s1.dec)**3 )
mask = err <= epsilon.to(u.radian).value
if (dec==45):
    t = lam/(2*np.cos(s1.dec))
else:
    t = (-np.cos(s1.dec) + np.sqrt(np.cos(s1.dec)**2 + lam**2 * np.cos(2*s1.dec)))/(lam * np.cos(2 * s1.dec))
    t[mask] = lam[mask]/(2*np.cos(s1.dec))
adj_ra2 = 2 * np.arctan2(t.to(u.dimensionless_unscaled), 1)
err2 = (-adj_ra2 +s2.ra).to(u.arcsec)

# fig, ax = plt.subplots()
plt.scatter(s2.ra.degree, err1, marker = '.', label="cos(dec) correction")
# ax2 = ax.twinx()
plt.scatter(s2.ra.degree, err2, marker = '.', label="true correction")
plt.xlabel("Offset in RA (in degree)", fontsize=15)
plt.ylabel("Error in RA correction (arcsec)", fontsize=15) 
plt.title(f"Fields at declination = {s1.dec.degree}", fontsize=15)
plt.legend(fontsize=15)
plt.tight_layout()

`

Note: I have to put an artificial cut to avoid numerical overflow errors for very small numbers

[dlakaplan]

Is the error really 10" for a 5-deg field at Dec=45?

[AkashA98]

Is the error really 10" for a 5-deg field at Dec=45?

Yes (edge-to-edge 5 degrees), redid my math, rechecked my code, and it seems correct, and it is roughly consistent with what I was seeing earlier (~ degree error over 30-40 degree offset in RA --- I recover this when I extend it to such offsets)

[dlakaplan]

10" would be much larger an error than we want. But I'm not sure if that's relevant, since I think this is just correcting the image center and all subsequent coordinate transformations are handled in the tangent pane.

[dlakaplan]

So I am pretty sure this is not an issue. Here are some reasons/tests:

• Take an image with the associated selavy catalog. Offset the image header by (Dra,Ddec) using the previous method (Dra/cos(Dec) etc). Offset the catalog properly in astropy. Make a region file from the new catalog and plot on the updated image. Do the sources line up? Any biases as a function of Dec? I am pretty sure Andrew & I did this test when developing code.

• Find an image with a significant shift wrt the catalog. Correct the header, remeasure the sources, and compare against the catalog again. Has the shift converged to 0 or is there any Dec-dependent error? Again, I think Andrew tested this.

Also, I may have misinterpreted your plot. What is relevant is the error in RA for a ~10" shift, but over ~+/-2.5 degrees in Dec.

[AkashA98]

So I am pretty sure this is not an issue. Here are some reasons/tests:

• Take an image with the associated selavy catalog. Offset the image header by (Dra,Ddec) using the previous method (Dra/cos(Dec) etc). Offset the catalog properly in astropy. Make a region file from the new catalog and plot on the updated image. Do the sources line up? Any biases as a function of Dec? I am pretty sure Andrew & I did this test when developing code.
• Find an image with a significant shift wrt the catalog. Correct the header, remeasure the sources, and compare against the catalog again. Has the shift converged to 0 or is there any Dec-dependent error? Again, I think Andrew tested this.

Also, I may have misinterpreted your plot. What is relevant is the error in RA for a ~10" shift, but over ~+/-2.5 degrees in Dec.

Ok, yes, agree that for the offsets that we are dealing with, we might not need to complicate things (as I am doing), I'll revert the code to its older state. Thanks...

[ddobie]

@dlakaplan is correct - the maximum offset we're correcting for is <10" (typically something like 2-3") so based on your plots the error is negligible.

[AkashA98]

Yes, thanks @dlakaplan @ddobie . I changed the code and also fixed the error propagation in my PR. Please take a look at it and let me know if any further changes are needed before merging. Thanks...