Looking at afw.image.TanWcs.cc, we see that the TanSip decoding is ignoring the first terms of the A and B polynomials (constant and linear terms). At first look, this is logical because these terms are supposed to be in the CD matrix and the CRPIX, but the simultaneous astrometry is fitting general polynomials without assuming the SIP representation. Converting from these polynomials to a standard SIP is not trivial and lead to numerical instabilities unless we use the linear and constant terms.
Apparently other SIP decoding implementations like : wcstools-3.9.1/libwcs/wcssubs-3.9.1/distort.c are not ignoring the first SIP terms.
Looking at the document describing the SIP convention it is not specified that these terms should be explicitly ignored and they should be taken into account for the reverse transformation (AP / BP).
So the question is : “should we modify TanWcs.cc to take into account the linear and constant terms of the TanSip if they are non zero” ? The change is trivial : just changing the lower bound of the loop index.
If other SIP implementations read those terms and the convention description is vague, it sounds like we probably should too.
I’m a bit concerned this would make the SIP representation non-unique (two exactly equivalent transforms could be written in different ways), but I don’t know how much of a problem that is (or even whether it’s already possible).
Look at equation 1: “Let f(u,v) and g(u,v) be the quadratic and higher- order terms of the distortion polynomial”
So SIP requires that you adjust the polynomials so that ignoring A/B reduces to the TAN projection.
Is this crazy? Probably. Is SIP a bad standard, using raw polynomials? Yes. But it’s one of the “standards” that people use. I am strongly in favour of overhauling our astrometric distortion approach (but that runs into questions of AST that we haven’t resolved).
Once the polynomial transformation has been fitted, in can be spllt in several ways between the standard
linear part and the SIP correction. One way is indeed to force the SIP low-order terms to zero.
It means that we have to find a point in pixel space where P(X)=0. This point will become the reference
pixel (CRPIX), and the linear expansion of the polynomial at this point will provide the CD terms, and the remainder provides the SIP A and B (free by construction of constant and linear terms). The problem is that this point where P(X) is zero is not necessarily unique, and can be far away from the bounding bow where the polynomial makes sense. Indeed, this distance depends on the choice of the reference point. The obvious alternative is to split the fitted polynomial by expanding around the center of the chip, because doing so, the SIP corrections are small.
In the end, the whole thing boils down to how we chose the reference points of WCSs. If we want
to be able to take the same reference point for a whole exposure, we will run into problems by Taylor
expanding our corrections at the wrong point. Polynomial corrections are not supposed to make sense
far beyond the domain they were fitted on.
In any case, one change has to be applied to afw: the constant and linear terms of the A and B polynomials are read from the headers, provided by accessors, but ignored when it comes to transforming. This is inconsistent.
And if the plan is to not allow constant and linear terms in A and B sip polynomials, then
meas-simastrom has to change significantly: we have to give up the scheme where
a whole exposure uses the same reference point.
I think the plan for now is to solve globally, then derive an approximate TAN-SIP for each small region (with different CRPIX points).
I totally agree that a sane internal focal plane description is a very high priority, keeping TAN-SIP and TAN-TPV as export formats for external codes.
I had a chat about this issue with @boutigny and @PierreAstier at the DESC meeting earlier this week, and I came away convinced that we need to support at least reading the constant and linear terms in afw, even it’s not permitted by the convention. Other popular codes that parse SIP (at least wcstools and DS9) do read them, and I think from a usability standpoint it’s better to be permissive than proscriptive; it’s not like there’s any confusion about how those terms should be interpreted.
More importantly, reading them will let us move ahead with more interesting aspects of meas_simastrom while we work on an afw Wcs class that can handle more interesting (i.e. multi-element AST-style) transformations natively. Once we have that available, meas_simastrom can be upgraded to fit more interesting models and we can relegate the SIP representation to just an approximate version we persist with the true transformation for convenience. At that point I’d be a lot more sanguine about automatically transforming illegal SIPs into legal ones when meas_simastrom writes them.