Description
Hello, I am a julia rookie, so I lost a fair time to discover that the inconsistent results I got from a FITS file produced by a colleague were due to using the WCS routines (the problem probably can lie in the header of the file but I am not expert in this)
The fits file contains a 4096 x 2048 array to be plotted. The projection is such that RA grows "to left" in pix space but it is wrapped around about 288 RA deg for projection purposes (so the origin RA=0 is not centered at 2048.5). DEC range is the standard [-90,90]. Exact WCSTransform is accluded below.
Both going from world coords to pixels or back has problems.
I acclude results from one simple tests, comparing them with an approximate "by hand" transform from RA to pix that the colleague gave me (the latter is off by 0.5 deg but this is unimportant here) and the correct results by using python astropy6.0.0 (the outputs in this case have RA <0, so need to wrap around them by adding 360)
Oddily, the Julia WCS routine yields mostly null outputs.
I suspect that the WCS from the file might miss some info that can avoid the incorrect behaviour but I have no idea about and I puzzle why astropy instead gives correct results.
---- failing WCSTransform ----
after reading from the file, the full header is
SIMPLE = T
BITPIX = -64
NAXIS = 2
NAXIS1 = 4096
NAXIS2 = 2048
EXTEND = T / FITS dataset may contain extensions
COMMENT FITS (Flexible Image Transport System) format is defined in 'Astronom
COMMENT and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H
CTYPE1 = 'RA---CAR'
CTYPE2 = 'DEC--CAR'
CRVAL1 = -252.0
CRVAL2 = 0.0
CRPIX1 = 2048.5 / Reference Pixel in X
CRPIX2 = 1024.5 / Reference Pixel in Y
CDELT1 = -0.087890625
CDELT2 = 0.087890625
CROTA2 = 0.0
once the above header is read as a single string (please stress this point in the info, it took me a while to find out why the header if read without reading it with the String keyword was failing)
I get:
-------------------- (culprit?) WCSTransform
wcshere=WCSTransform(naxis=2,cdelt=[-0.066667, 0.066667],crval=[0.0, -90.0],crpix=[-234.75, 8.3393])
------------------ test code - (copied form notebook, output below)
WCS.wcslib_version()
fitsfilename="XXX.fits"
fitsfile = FITS(fitsfilename)
fileheaderstring = FITSIO.read_header(fitsfile[1],String)
println(fileheaderstring, "\n")
wcshere=WCS.from_header(fileheaderstring)[1]
println(wcshere, "\n")
# test where origin RA=0 and DEC=0 goes in pixel space
radecorigin=zeros(2,1)
# the transform for origin seems to be OK
pixorigin=WCS.world_to_pix(wcshere, radecorigin )
backradecorigin =WCS.pix_to_world(wcshere, pixorigin )
println("\n pixorigin ",pixorigin," backradecorigin",backradecorigin )
# test values
vra = [0.0, 90.0, 180.0, 270, 359.0, 360.0, 0.0, 0.0, 180.0, 359.0, 360.0, 360.0]
vdec = [0.0, 0.0, 0.0,0.0, 0.0, 0.0, 89.0, 90.0, 90.0, 90.0, 90.0,0.0]
radecmatorig=[ vra ;; vdec]
println(typeof(radecmatorig))
println("radecmatorig \n",radecmatorig)
display(radecmatorig)
sz=size(radecmatorig)
npoi=sz[1]
println(" sz ",sz, "\n npoints ",npoi, )
pixmat = WCS.world_to_pix(wcshere, radecmatorig )
println(" pixmat \n",pixmat)
backradec=WCS.pix_to_world(wcshere, pixmat)
# simple "by hand" transform
naxis1=4096
naxis2=2048
deg2pix=float(naxis1)/360.0
handpixmat=similar(pixmat)
for k= 1:npoi
ra=radecmatorig[k,1]
dec=radecmatorig[k,2]
if ra > 288
x = naxis1 - (ra - 288) * deg2pix
else
x = (288 - ra) * deg2pix
end
y = dec * deg2pix + naxis2 / 2
handpixmat[k,1]=x
handpixmat[k,2]=y
end
println("\n")
println("\n radecmatorig ")
display(radecmatorig)
println("\n pixmat ")
display(pixmat)
println("\n handpixmat ")
display(handpixmat)
println(typeof(radecmatorig))
println("radecmatorig \n",radecmatorig)
display(radecmatorig)
sz=size(radecmatorig)
npoi=sz[1]
println(" sz ",sz, " npoi \n",npoi, )
---- julia output ------
SIMPLE = T BITPIX = -64 NAXIS = 2 NAXIS1 = 4096 NAXIS2 = 2048 EXTEND = T / FITS dataset may contain extensions COMMENT FITS (Flexible Image Transport System) format is defined in 'AstronomyCOMMENT and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H CTYPE1 = 'RA---CAR' CTYPE2 = 'DEC--CAR' CRVAL1 = -252. CRVAL2 = 0. CRPIX1 = 2048.5 / Reference Pixel in X CRPIX2 = 1024.5 / Reference Pixel in Y CDELT1 = -0.087890625 CDELT2 = 0.087890625 CROTA2 = 0.0 END
WCSTransform(naxis=2,cdelt=[-0.087890625, 0.087890625],crval=[-252.0, 0.0],crpix=[2048.5, 1024.5])
pixorigin [3277.3; 1024.5;;] backradecorigin[-1.4210854715202004e-14; 0.0;;]
radecmatorig
12×2 Matrix{Float64}:
0.0 0.0
90.0 0.0
180.0 0.0
270.0 0.0
359.0 0.0
360.0 0.0
0.0 89.0
0.0 90.0
180.0 90.0
359.0 90.0
360.0 90.0
360.0 0.0
sz (12, 2)
npoi 12
## BAD OUTPUT
pixmat
12×2 Matrix{Float64}:
3277.3 3277.3
2048.5 1024.5
2.36776e-314 2.36776e-314
0.0 0.0
2.36775e-314 2.36775e-314
2.36776e-314 2.36776e-314
0.0 0.0
2.36775e-314 2.36775e-314
2.35767e-314 2.36776e-314
0.0 0.0
2.36775e-314 5.31645e-314
2.36776e-314 5.31645e-314
radecmatorig
12×2 Matrix{Float64}:
0.0 0.0
90.0 0.0
180.0 0.0
270.0 0.0
359.0 0.0
360.0 0.0
0.0 89.0
0.0 90.0
180.0 90.0
359.0 90.0
360.0 90.0
360.0 0.0
## BAD BAD OUTPUT
pixmat
12×2 Matrix{Float64}:
3277.3 3277.3
2048.5 1024.5
2.9969e-314 2.36776e-314
0.0 0.0
2.36775e-314 2.36775e-314
2.36776e-314 2.36776e-314
0.0 0.0
2.36775e-314 2.36775e-314
2.35767e-314 2.36776e-314
0.0 0.0
2.36775e-314 5.31675e-314
2.36776e-314 5.31675e-314
handpixmat (always off by 0.5 degs)
12×2 Matrix{Float64}:
3276.8 1024.0
2252.8 1024.0
1228.8 1024.0
204.8 1024.0
3288.18 1024.0
3276.8 1024.0
3276.8 2036.62
3276.8 2048.0
1228.8 2048.0
3288.18 2048.0
3276.8 2048.0
3276.8 1024.0
backradec ### BAD OUTPUT ###
12×2 Matrix{Float64}:
-1.42109e-14 -1.42109e-14
90.0 0.0
0.0 0.0
0.0 0.0
0.0 0.0
0.0 0.0
0.0 0.0
0.0 0.0
0.0 0.0
0.0 0.0
0.0 0.0
0.0 0.0
------ same conversion done successfully in python
....... imports, reads etc....
print(mradec)
pixmat=w.wcs_world2pix(mradec,1)
print(pixmat)
backmat=w.wcs_pix2world(pixmat,1)
print,backmat
# wrap around negative RA coords
negVals = np.where(backmat[:,0] < 0.)
print(negVals)
print(backmat)
backmat[negVals,0] += 360.0
print(mradec)
print(backmat)
---- output # Correct
# ra dec array
[[ 0. 0.]
[ 90. 0.]
[180. 0.]
[270. 0.]
[359. 0.]
[360. 0.]
[ 0. 89.]
[ 0. 90.]
[180. 90.]
[359. 90.]
[360. 90.]
[360. 0.]]
# pixel array
[[3277.3 1024.5 ]
[2253.3 1024.5 ]
[1229.3 1024.5 ]
[ 205.3 1024.5 ]
[3288.67777778 1024.5 ]
[3277.3 1024.5 ]
[3277.3 2037.12222222]
[3277.3 2048.5 ]
[1229.3 2048.5 ]
[3288.67777778 2048.5 ]
[3277.3 2048.5 ]
[3277.3 1024.5 ]]
# back transformed ra dec array
array([[-1.42108547e-14, 0.00000000e+00],
[-2.70000000e+02, 0.00000000e+00],
[-1.80000000e+02, 0.00000000e+00],
[-9.00000000e+01, 0.00000000e+00],
[-1.00000000e+00, 0.00000000e+00],
[-1.42108547e-14, 0.00000000e+00],
[-1.42108547e-14, 8.90000000e+01],
[-1.42108547e-14, 9.00000000e+01],
[-1.80000000e+02, 9.00000000e+01],
[-1.00000000e+00, 9.00000000e+01],
[-1.42108547e-14, 9.00000000e+01],
[-1.42108547e-14, 0.00000000e+00]]))
### negative RA wrap around done here
# output array before back transform and wrap
(array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]),)
[[-1.42108547e-14 0.00000000e+00]
[-2.70000000e+02 0.00000000e+00]
[-1.80000000e+02 0.00000000e+00]
[-9.00000000e+01 0.00000000e+00]
[-1.00000000e+00 0.00000000e+00]
[-1.42108547e-14 0.00000000e+00]
[-1.42108547e-14 8.90000000e+01]
[-1.42108547e-14 9.00000000e+01]
[-1.80000000e+02 9.00000000e+01]
[-1.00000000e+00 9.00000000e+01]
[-1.42108547e-14 9.00000000e+01]
[-1.42108547e-14 0.00000000e+00]]
# orignal input array
[[ 0. 0.]
[ 90. 0.]
[180. 0.]
[270. 0.]
[359. 0.]
[360. 0.]
[ 0. 89.]
[ 0. 90.]
[180. 90.]
[359. 90.]
[360. 90.]
[360. 0.]]
# backtransformed array
[[360. 0.]
[ 90. 0.]
[180. 0.]
[270. 0.]
[359. 0.]
[360. 0.]
[360. 89.]
[360. 90.]
[180. 90.]
[359. 90.]
[360. 90.]
[360. 0.]]