These are chat archives for astropy/astropy

May 2016
May 21 2016 10:54
My WCS quest continues. After Calabretta helpfully showed me how to use LONPOLE to set the obliqueness independently of the reference point, I have a nicely working example in C. However, the equivalent version in python acts as if lonpole is simply being ignored
These should produce the same output
In particular, they should show declination (the last column) changing in steps of 10 for each row, while ra stays mostly the same
May 21 2016 11:10
Yet the first few lines of output from the C version is:
0 0 0 1.000 1.000 360.000 0.000
3 0 0 4.000 1.000 360.000 30.000
4 0 0 5.000 1.000 360.000 40.000
1 0 0 2.000 1.000 360.000 10.000
2 0 0 3.000 1.000 360.000 20.000
5 0 0 6.000 1.000 360.000 50.000
which is correct
while the corresponding lines from the python version are:
0 1.000 1.000 360.000 0.000
3 4.000 1.000 330.000 0.000
2 3.000 1.000 340.000 0.000
4 5.000 1.000 320.000 0.000
1 2.000 1.000 350.000 0.000
5 6.000 1.000 310.000 0.000
May 21 2016 11:24
I can reproduce the python version's behavior by setting latpole = -90 in the C version
And inspection of the python wcs object shows that indeed latpole is -90 there
When I try to change latpole, it seems to switch to either +90 or -90, no matter what value I actually set
(at least while lonpole = 90)
May 21 2016 12:12
This is weird. latpole does get properly set in the wrapper. But then later, when it needs to be used, the value has somehow changed
May 21 2016 12:58
It seems to be overriden no matter what I do
It happens in wcs.set(), and in wcs.bounds_check() and in wcs.p2s
And it always happens in wcslib code
But when use wcslib directly in C, it doesn't happen
so something must be different
May 21 2016 13:25
Wow, this difference is pretty obscure!
It decides whether to override LATPOLE or not based on the latpre variable
This one is 2 for the working C case, and one for the non-working python case
The check that determines which of these will happen is on line 285 in cel.c in wcslib: "if (z == 0.0)"
In the C case, z is 0. In the python case, z is 6.123233995736766e-17
So this looks like some sort of floating point precision issue
But wait! There are larger differences before this point
z gets is value from x, which has the same value. And x = cthe0*cphip.
cthe0 = 1, while cphip has the value 0 for C and 6.123233995736766e-17 for python
cphip is computed by sincosd(phip - cel->phi0, &sphip, &cphip)
May 21 2016 13:30
And interestingly, it looks like the unit of the output is different for the two versions
Actually, strike the last sentence
Everything seems the same until the point where cphip becomes 0 in C and a tiny nonzero value in python
May 21 2016 13:53
This difference seems to come from which version of sincosd is being called
In the C version, the function defined in wcstrig.c is called. This one has checks for poles, etc. In the python version, a simple macro is used instead. This macro does not have those checks
Which to use is defined by the preprocessor symbol WCSTRIG_MACRO, which I assume is set for the python version but not the C version
May 21 2016 14:00
It looks like astropy defines this macro. Trying a compile without it now
Looks like wcstrig.c isn't even compiled in astropy
May 21 2016 14:05
Ok, after adding wcstrig.c to the files to be compiled, and disabling WCSTRIG_MACRO, the python version finally works
May 21 2016 14:24
Bug report sent to Calabretta
The real bug is in wcslib, but astropy can easily work around it by changing two lines in astropy/wcs/
Erik Tollerud
May 21 2016 16:42
@amaurea - sounds reasonable - I'd encourage you to submit a PR to make that change yourself! (If you're not familiar with the process, a good start is
@MSeifert04 - when you say "position" for do you mean which it's in, or where in the same file. The former does matter (and conceivably could be different for RTD vs. local... although it "shouldn't"), but the latter should always be the same, because is parsed as just a regular python module (on RTD or locally).
May 21 2016 17:18
@eteq: Can I update a pull request on github?
I need to add tests to the other one, I guess
May 21 2016 17:52
Bah, I'll need to make a new pull request
I didn't make a branch for the previous changes, and I've since added another feature
Matt Craig
May 21 2016 18:10
@amaurea yes, you can update a pull request by pushing additional commits to it...But I agree a new PR on a branch is much, much easier to manage!