[AstroPy] WCS query

Stephen Walton stephen.walton at csun.edu
Wed Jun 15 12:02:29 CDT 2005


Perry Greenfield wrote:

> Hi Steve,
>
> I have a feeling this is related to your matplotlib query.

Yes.

>
> On Jun 10, 2005, at 7:57 PM, Stephen Walton wrote:
>
>> However, if I read the WCS with PyFITS and try to compute a WCS array 
>> in Python, things go wrong.  Here's what I did:  I created a 2x2 cd 
>> array in Python where cd[k-1,m-1] was set to CD_km from the FITS 
>> header; I also read CRPIX1 and CRPIX2 from the FITS header into 
>> crpix1 and crpix2 variables.  Two arrays, x and y, were produced, 
>> where x[k,m]=k and y[k,m]=m, k,m=0,1023.  Finally, I created xx and 
>> yy as follows:
>>
> The above suggests you were associating x coordinates with the first 
> index and y with the second. Of course with numarray (and Numeric) 
> that's wrong. Is this the crux of the problem?

As it turns out, yes, and the code below looks correct and seems to 
produce the right result:

x,y=meshgrid(range(1024),range(1024))   # x[i,j] == j, y[i,j] == i
xx=cd[0,0]*(x-crpix1+1)+cd[0,1]*(y-crpix2+1)
yy=cd[1,0]*(x-crpix1+1)+cd[1,1]*(y-crpix2+1)

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.

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.




More information about the AstroPy mailing list