[AstroPy] WCS query

Perry Greenfield perry at stsci.edu
Wed Jun 15 16:02:54 EDT 2005

On Jun 15, 2005, at 1:02 PM, Stephen Walton wrote:

> But I still find it somewhat awkward, and users here are going to find 
> it very confusing, that what we old Fortran hands think of as FITS 
> pixel (i,j) winds up displayed in matplotlib at Cartesian coordinates 
> (i-1,j-1) even though that pixel would have to be addressed as 
> data[j-1,i-1] in software.  This is because John Hunter, not 
> unreasonably, adopted the MATLAB convention for imshow(), contour(), 
> and their relatives that the first array coordinate is vertical and 
> the second horizontal.  Effectively, then, imshow() 'undoes' the 
> transpose done by PyFITS so that the displayed image actually comes 
> out in the same orientation on the screen as it does in IRAF (if 
> origin='lower' is used in imshow, as I do) even though the actual 
> array has to be addressed with the indices in reverse order from 
> Fortran.

This whole area hinges on the definitions of what the conventions are 
(which can be confusing). I'm not sure I'd classify it as John Hunter's 
convention but more properly a Numeric/numarray convention. That's 
because the data as read from the FITS file retain the same order they 
did in the FITS file (PyFITS doesn't reorder the data). Since the last 
index for Numeric/numarray represents the adjacent data values, that's 
where the issue really arises, not in matplotlib other than it is 
adopting the convention that adjacent data values are associated with x 
as far as images go (but having two different conventions for which way 
y should be displayed, i.e., up or down).

As I've mentioned, this is probably the biggest sore point astronomers 
are going to have with numarray (or Numeric). Different approaches 
could have been used, but all have their drawbacks:

1) numarray should have adopted the other convention. It could have, 
but then it would have been incompatible with Numeric in this respect 
and that would have been bad. [Numeric chose (I believe) to do it that 
way since it means that arr[i,j] == arr[i][j] instead of 
arr[i,j]==arr[j][i] and that was thought important enough to be 
consistent for]

2) PyFITS could actually reorder the data on reading. I suppose, but 
this can be an expensive operation for large datasets and it renders 
memory mapping useless. Even after doing so, the presumption that x 
corresponds to adjacent data values is no longer true (so if algorithms 
are written that presume that, large inefficiencies will result)

3) PyFITS could read the data in the normal order and fiddle with the 
array stride attributes to have the same effect. This would allow 
continued use of memory mapping. Although now the first index  
represents adjacent data, all other Numeric/numarray libraries that 
presume the last index usually represents adjacent data will be wrong.

Any of these choices results in some sort of mismatch. We concluded 
that, although initially painful, the best thing was to remain as 
consistent with the existing language and array package semantics 
rather than introduce tricks to make it look like another language's 
semantics. The former is much more painful up front, the latter 
introduces long-term problems that never go away. Perhaps that's a bad 
marketing strategy. It's a little like the 0-based vs 1-based indexing 
issue. One can try to set things up so that you can use 1-based 
indexing in a 0-based language, but in the end it's just more trouble 
than learning to use it the way the language works (in my opinion 
anyway). But I'll agree that this is probably the single biggest 
irritation that astronomers will face.

> Incidentally, MATLAB's provided fitsread routine transposes FITS data 
> on input, so in MATLAB data(i,j) refers to the same pixel as in 
> Fortran (and image displays *look* transposed).
> I suppose long-time CFITSIO and PyFITS users are just used to all of 
> this.
That depends. There was never consensus even in our group about which 
approach was best (though there did seem to be an age correlation, but 
maybe I imagined that).


More information about the AstroPy mailing list