Cartesian projection issue in a FITS image through PyFITS / AstroPy

Question

I've looked and looked for a solution to this problem and am turning up nothing.

I'm generating rectangular FITS images through matplotlib and subsequently applying WCS coordinates to them using AstroPy (or PyFITS). My images are in galactic latitude and longitude, so the header keywords appropriate for my maps should be GLON-CAR and GLAT-CAR (for Cartesian projection). I've looked at other maps that use this same map projection in SAO DS9 and the coordinates work great... the grid is perfectly orthogonal as it should be. The FITS standard projections can be found here.

But when I generate my maps, the coordinates are not at all Cartesian. Here's a side-by-side comparison of my map (left) and another reference map of roughly the same region (right). Both are listed GLON-CAR and GLAT-CAR in the FITS header, but mine is screwy when looked at in SAO DS9 (note that the coordinate grid is something SAO DS9 generates based on the data in the FITS header, or at least stored somewhere in the FITS file):

This is problematic, because the coordinate-assigning algorithm will assign incorrect coordinates to each pixel if the projection is wrong.

Has anyone encountered this, or know what could be the problem?

I've tried applying other projections (just to see how they perform in SAO DS9) and they come out fine... but my Cartesian and Mercator projections do not come out with the orthogonal grid like they should.

I can't believe this would be a bug in AstroPy, but I can't find any other cause... unless my arguments in the header are incorrectly formatted, but I still don't see how that could cause the problem I'm experiencing. Or would you recommend using something else? (I've looked at matplotlib basemap but have had some trouble getting that to work on my computer).

My header code is below:

 from __future__ import division
 import numpy as np
 from astropy.io import fits as pyfits # or use 'import pyfits, same thing'

 #(lots of code in between: defining variables and simple calculations...
 #probably not relevant)

 header['BSCALE'] = (1.00000, 'REAL = TAPE*BSCALE + BZERO')
 header['BZERO'] = (0.0)
 header['BUNIT'] = ('mag ', 'UNIT OF INTENSITY')
 header['BLANK'] = (-100.00, 'BLANK VALUE')
 header['CRVAL1'] = (glon_center, 'REF VALUE POINT DEGR')   #FIRST COORDINATE OF THE CENTER
 header['CRPIX1'] = (center_x+0.5, 'REF POINT PIXEL LOCATION') ## REFERENCE X PIXEL
 header['CTYPE1'] = ('GLON-CAR', 'COORD TYPE : VALUE IS DEGR')
 header['CDELT1'] = (-glon_length/x_length, 'COORD VALUE INCREMENT WITH COUNT DGR')   ### degrees per pixel
 header['CROTA1'] = (0, 'CCW ROTATION in DGR')            
 header['CRVAL2'] = (glat_center, 'REF VALUE POINT DEGR') #Y COORDINATE OF THE CENTER
 header['CRPIX2'] = (center_y+0.5, 'REF POINT PIXEL LOCATION') #Y REFERENCE PIXEL 
 header['CTYPE2'] = ('GLAT-CAR', 'COORD TYPE: VALUE IS DEGR')   # WAS CAR OR TAN
 header['CDELT2'] = (glat_length/y_length, 'COORD VALUE INCREMENT WITH COUNT DGR') #degrees per pixel  
 header['CROTA2'] = (rotation, 'CCW ROTATION IN DEGR')                               #NEGATIVE ROTATES CCW around origin (bottom left). 
 header['DATAMIN'] = (data_min, 'Minimum data value in the file') 
 header['DATAMAX'] = (data_max, 'Maximum data value in the file')
 header['TELESCOP'] = ("Produced from 2MASS")

 pyfits.update(filename, map_data, header)

Thanks for any help you can provide.

Answer 1

In the modern definition of the -CAR projection (from Calabretta et al.), GLON-CAR/GLAT-CAR projection only produces a rectilinear grid if CRVAL2 is set to zero. If CRVAL2 is not zero, then the grid is curved (this should have nothing to do with Astropy). You can try and fix this by adjusting CRVAL2 and CRPIX2 so that CRVAL2 is zero. Does this help?

Just to clarify what I mean, try, after your code above, and before writing out the file:

header['CRPIX2'] -= header['CRVAL2'] / header['CDELT2']
header['CRVAL2'] = 0.

Any luck?

If you look at the header for the 'reference' file you looked at, you'll see that CRVAL2 is zero there. Just to be clear, there's nothing wrong with CRVAL2 being non-zero, but the grid is then no longer rectilinear.