search

JuliaAstro / WCS.jl

Astronomical World Coordinate Systems library for Julia
https://juliaastro.github.io/WCS.jl/latest/
MIT License
16 stars 12 forks source link

Failure of an unusual WCS which, however, works on python #47

Open Illuminista opened 8 months ago

Illuminista commented 8 months ago

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.]]
astrozot commented 8 months ago

The issue is related to the use of 12 x 2 matrix for the list of coordinates (radecmatorigabove) instead of a 2 x 12 matrix. This creates the wrong result in the world_to_pix transformation.

Essentially, this kind of mistakes happen because the library is quite a thin wrapper to the C WCS library, with very little checks for input parameters.