These are chat archives for astropy/astropy

18th
May 2016
amaurea
@amaurea
May 18 2016 15:40
Ok, trying again in this channel [yay for proprietary, centralized alternatives to irc]
Are there any WCS experts here? I'm trying to define a cylindrical coordinate system that follows a line of longitude.
I've got it almost working by setting CRVAL2 = 90
The problem is that that also sets the reference point to 90°
And wcslib has this nasty habit (bug?) of refusing to transform points that are further than 180° away from the reference point. I can usually work around this by moving the reference point such that my entire image falls within the allowed range. But how do I do that in this case?
I've messed around with LATPOL, LONPOL, PV1_X and PV0_1, but the standard is confusing and I haven't gotten anything to work
Erik Tollerud
@eteq
May 18 2016 16:36
@amaurea - what projection are you trying to do this with?
@amaurea (also, there’s an IRC bridge to gitter that you can use if you prefer IRC: https://irc.gitter.im/)
amaurea
@amaurea
May 18 2016 16:39
eteq: I'm trying it with CAR
I guess there are two issues here, really. The first is the NaNs you get when you wrap around the sky. This happens even in perfectly normal coordinates
The second one is then the ability to choose the reference point without changing the projection
Fixing any of these would work, but the first one would be the most convenient, as I've encountered it other places too
I've looked a bit more into this now, and interestingly the blame seems to be shared between wcslib and astropy
amaurea
@amaurea
May 18 2016 16:44
Wcslib itself does actually properly wrap around the sky multiple times (at least in my simple test). But after it gets too far away from the reference point, it sets the status return value to 8
astropy then sees this, and overwrites the (actually correct) output from wcslib with NaN
Erik Tollerud
@eteq
May 18 2016 16:45
One possible “quick fix” to that: astropy.coordinates.Angle(angleval, angleuint).wrap_at(180*u.deg)
that’s a bit hacky, but should do the trick, right?
amaurea
@amaurea
May 18 2016 16:46
The problem is in pix2sky, though
Erik Tollerud
@eteq
May 18 2016 16:46
(that is, you do that before passing it into wcslib)
amaurea
@amaurea
May 18 2016 16:46
So you would have to wrap the pixels instead
Erik Tollerud
@eteq
May 18 2016 16:46
ohhh
wait, but how do you get much greater than 180 degrees with that mapping? That would imply you’ve got a detector that wraps around the sky twice...
amaurea
@amaurea
May 18 2016 16:47
I'll give you a code example
(but what you say would only apply if the reference point is in the center of your field of view. It might not be, and you may not be allowed to place it there by the projection)
Erik Tollerud
@eteq
May 18 2016 16:49
I admit I don’t really understand the starting point (i.e. what you meant by “trying to define a cylindrical coordinate system that follows a line of longitude”), but maybe your example will help with that
amaurea
@amaurea
May 18 2016 16:50
Yes, I will give you two examples, first a plain one with normal car, then a rotated one
But I stupidly deleted my example by mistake, so I'm re-typing it
amaurea
@amaurea
May 18 2016 16:56
The reference point is at ra=0,dec=0, so pixels that would result in ra>180 result in nan
But in this case one could just move the reference point along the equator to work around this. Moving the reference point along the equator does not change the projection
eteq: If you run it, you should see sensible output for the first few points, and then nan for the rest
Here's the example with a rotated CAR system: http://folk.uio.no/sigurdkn/wcswrap/wrap2.py
Here the "equator" goes along the ra=0 line of longitude
This system is good for representing the set of images a satellite in polar orbit around would take
Again things work fine while one is near the reference point
amaurea
@amaurea
May 18 2016 17:02
But this time we don't have the freedom of moving it
If we move it away from dec=90, then we switch to a different projection that no longer follows a line of longitude
I also have C example programs that do exactly the same thing
But these interestingly show that the wcslib output is correct, except that it needlessly indicates that the output is invalid
Here is a C program corresponding to the second example: http://folk.uio.no/sigurdkn/wcswrap/wrap2.c
amaurea
@amaurea
May 18 2016 17:41
If wcslib keeps doing this, perhaps it would be possible to add an option to wcs_pix2world and wcs_world2pix to make them ignore the wcslib return status, and treat all the output as valid?
amaurea
@amaurea
May 18 2016 17:55
In WCSlib, the particular error code BAD_PIX is special-cased to let the calculations continue
The error code is set in prjbchk in prj.c in wcslib. It has explicit checks about -180 < phi < 180. I guess they're there for a reason, but I would be tempted to just remove them
amaurea
@amaurea
May 18 2016 18:05
(the reason why they're there is probably to make world2pix(pixworld(x)) == x at all times. I'm not sure that's worth it, though
@eteq: Oh, I guess one has to prefix by @ here when referencing someone
Erik Tollerud
@eteq
May 18 2016 18:13
(sorry, had something come up, @amaurea - will try to have a closer look in a bit)
amaurea
@amaurea
May 18 2016 18:17
The wraparound issue discussed on page 36 in part 2 of the standard: http://www.atnf.csiro.au/people/mcalabre/WCS/ccs.pdf
As far as I can see, it does not say that going more than 180° away from the reference point is illegal
amaurea
@amaurea
May 18 2016 18:29
back in 30 min
amaurea
@amaurea
May 18 2016 18:54
back
Erik Tollerud
@eteq
May 18 2016 19:06
Can you point to where this is mentioned in the WCS II paper?
Also, @amaurea, note that you can make code snippets in gitter with triple-backticks like this:
import numpy as np
from astropy.wcs import WCS
wcs = WCS(naxis=2)
wcs.wcs.ctype = ["RA---CAR","DEC--CAR"]
wcs.wcs.cdelt = [10,10]
wcs.wcs.crval = [0,0]
wcs.wcs.crpix = [1,1]

for i in range(40):
    x, y = i, 0
    ra, dec = wcs.wcs_pix2world(x,y,0)
    print "%2d %8.3f %8.3f %8.3f %8.3f" % (i, x, y, ra, dec)
amaurea
@amaurea
May 18 2016 19:08
ok
eteq: It is mentioned on page 36, in example 3 of how to interpret a fits header
They have an example there of a wide image where the reference point is set such that the image straddles the ±180 degree cut
The page is spent discussing how this makes things more difficult, and how to move the reference point to avoid this
amaurea
@amaurea
May 18 2016 19:22
@eteq: Can you reproduce the behavior in my examples?
Erik Tollerud
@eteq
May 18 2016 19:22
gotcha. It seems to imply that this is not an intended use case, which is probably why it was given an error code… but I see your point taht it might be useful
amaurea
@amaurea
May 18 2016 19:23
Yes, especially for example 2. Though perhaps there are some obscure keywords that would make that work too
I've e-mailed the wcslib author to ask him about it
Erik Tollerud
@eteq
May 18 2016 19:23
(I still don’t understand why you want a reference point of dec=90, but I’m not sure that’s especially important for this situation)
amaurea
@amaurea
May 18 2016 19:23
Oh, I'll try to explain that again
I don't want a reference point at dec=90. I'm forced to put it there
In normal CAR, the equator of the coordinate system follows the line dec=0, right?
I want a rotated version of CAR where the equator of the coordinate system follows the line ra=0
We have an instrument that scans along such lines of longitude. CAR and other cylindrical projections are great for representing areas near great circles on the sky with little distortion, so they seem perfect for this. But the great circle in question is one that goes through the poles, not along dec=0
amaurea
@amaurea
May 18 2016 19:28
The only way of achieveing that I know of is to put the reference point at one of the poles
But if you know another way, that would be great
@eteq: Did that help?
Erik Tollerud
@eteq
May 18 2016 19:34
I think so… I’m just trying to decide if we can get by all this by finding a better projection for you
although now that I think about it, all projections will have the problem you’re talking about because you want to wrap multiple scans, right?
amaurea
@amaurea
May 18 2016 19:38
eteq: Multiple scans would be nice, but I don't need that. I don't need to cover more than 360° total
The rotated CAR can support a full 360°, but the pixels need to start at one of the poles to avoid wrapping in that case
I guess I could live with that. In the worst case I can remap all the pixels in the image
Perhaps I'm expecting too much flexibility from wcs.
Let's see what the wcslib author says
He's also one of the authors of the standard, after all
Erik Tollerud
@eteq
May 18 2016 19:49
yeah, I’d say the big question is whether he considers it “correct” or not that the NaNs are getting in there (i.e., should status code of 8 still be treated as “bad” even when it sort of makes sense here)
if not, then I think your proposal of adding a keyword at the python layer that lets the NaN-setting be supressed is a good solution
amaurea
@amaurea
May 18 2016 19:50
Right. The wcslib code itself treats that error condition different from the others
For the others, it stops as soon as the problem happens
But for that one, it keeps going, and simply flags the result in the end
That would put this issue in a grey zone, I guess
Erik Tollerud
@eteq
May 18 2016 19:51
right, and then in wcslib_wrap.c (line 1316 in current master), the “invalid” points are then set to nan
amaurea
@amaurea
May 18 2016 19:51
Yes. So one argument one could make is that astropy, too, should special case this particular error
Erik Tollerud
@eteq
May 18 2016 19:51
so it would be pretty trivial to change the status ==8 to status == 8 | keyword_set sort of thing
yes, exactly
amaurea
@amaurea
May 18 2016 19:52
Yes, exactly
:)
Erik Tollerud
@eteq
May 18 2016 19:52
:clap:
at any rate, lets see what mark says, but if he seems to think your use case should be supported, I’d encourage you to make a PR to astropy to do just that
amaurea
@amaurea
May 18 2016 19:53
Who is Mark? Calabretta?
Yes, I see, the M is for mark
Ok
He was just "M" in the papers, I wasn't sure
Erik Tollerud
@eteq
May 18 2016 19:57
oh, yeah, him