[bug #59459] Unclear WCS when both PC and CD exist in input, but conflict

[Mohammad Akhlaghi]

Update of bug #59459 (project gnuastro):

                  Status:                    None => Fixed                  
             Open/Closed:                    Open => Closed                 
                 Summary: Problem with SIP to TPV trasnformation => Unclear
WCS when both PC and CD exist in input, but conflict

    _______________________________________________________

Follow-up Comment #4:

With the help of Sachin, we think the problem has been fixed in Commit
7e74cba65d
<http://git.savannah.gnu.org/cgit/gnuastro.git/commit/?id=7e74cba65d>. But
please read below (it may be important for your project)...

In short, Sachin found out that this FITS file contains both the PC and CD
matrices of the WCS standard (which are different ways to map the pixels to
the celestial coordinates). And more importantly, Sachin recognized that they
produce different mappings! Usually a FITS file only has one of them, so this
was strange! When they exist together, they are usually just different
representations of the same mapping.

So I had a look at it (a small portion of it is shown below), and from the
position of the headers, I saw the cause of the problem: astrometry.net had
kept the original keywords but changed *some* of their names (except for the
'PCX_X' matric keywords!!!). Immediately after it, astrometry.net has just
added its finally derived WCS after it. Here is a part of the header:


                      / World Coordinate System (WCS)                         
 
COMMENT Original key: "WCSAXES"                                               
 
_CSAXES =                    2 / Number of coordinate axes                    
 
COMMENT Original key: "CRPIX1"                                                
 
_RPIX1  =              -1095.0 / Pixel coordinate of reference point
COMMENT Original key: "CRPIX2"                                                
 
_RPIX2  =               2914.0 / Pixel coordinate of reference point          
 
PC1_1   =   -0.999998140178468 / Coordinate transformation matrix element     
 
PC1_2   =   -0.001928636721656 / Coordinate transformation matrix element     
 
PC2_1   =   -0.001928636721656 / Coordinate transformation matrix element     
 
PC2_2   =    0.999998140178468 / Coordinate transformation matrix element     
 
COMMENT Original key: "CDELT1"                                                
 
_DELT1  =  6.22201157183109E-05 / [deg] Coordinate increment at reference
point 
                                                                              
 
                      / Versions and date                                     
 
DATE    = '2020-11-17T16:47:10' / file creation date (YYYY-MM-DDThh:mm:ss UT) 
 
CFITSIO = '3.47    '           / CFITSIO version.                             
 
WCSLIB  = '6.4     '           / WCSLIB version.                              
 
GSL     = '2.6     '           / GNU Scientific Library version.              
 
GNUASTRO= '0.13.50-4605'       / GNU Astronomy Utilities version.             
 
COMMIT  = '60b8c1a-dirty'      / Git's commit description in running dir.     
 
COMMENT Original key: "END"                                                   
 
COMMENT                                                                       
 
COMMENT --Start of Astrometry.net WCS solution--                              
 
COMMENT --Put in by the new-wcs program--                                     
 
COMMENT                                                                       
 
...
CUNIT1  = 'deg     ' / X pixel scale units                                    
 
CUNIT2  = 'deg     ' / Y pixel scale units                                    
 
CD1_1   =   -6.23207219845E-05 / Transformation matrix                        
 
CD1_2   =    9.60338675613E-07 / no comment                                   
 
CD2_1   =    3.66229401919E-07 / no comment                                   
 
CD2_2   =    6.28347154262E-05 / no comment                                   
 
IMAGEW  =                 2045 / Image width,  in pixels.                     
 
IMAGEH  =                 4605 / Image height, in pixels.                     
 
A_ORDER =                    2 / Polynomial order, axis 1                     
 
A_0_0   =                    0 / no comment                                   
 
A_0_1   =                    0 / no comment                                   
 
A_0_2   =   -2.97891115255E-06 / no comment                                   
 
A_1_0   =                    0 / no comment                                   
 
...


Because it had renamed the original CDELT keywords (which define the pixel
scale for the PC matrix), but had kept the PC matrix keywords, the two PC and
CD representations didn't match any more!!!

So I manually deleted the PC matrix elements with the command below, then
converted the distortion to TPV and the problem was fixed: I checked by
opening the TPV file with the SIP file without the PC matrix elements, in ds9.
Then I matched both of them by WCS and the images matched perfectly. See below
for the Python output that you originally posted.


cp in.fits in-copy.fits
astfits in-copy.fits --delete=PC1_1 --delete=PC1_2 --delete=PC2_1
--delete=PC2_2


Keeping old metadata (that is no longer relevant) is very bad practice, that
causes terrible situations like this. Unfortunately many astronomical software
follow this bad practice, causing headaches like this in later phases. But
anyway, let's get back to the problem.

So the problem wasn't actually caused the SIP --> TPV conversion! It was
caused by two different pixel <--> sky conversions in one file!!!

To avoid frustrations like this for future Gnuastro users, in Commit
7e74cba65d
<http://git.savannah.gnu.org/cgit/gnuastro.git/commit/?id=7e74cba65d>, I added
a check for such cases: if the two mappings exist in the file, and they are
different, Gnuastro will print a detailed warning and even suggest solutions.
Then it will use the CD matrix (as described in the commit message).

THIS IS IMPORTANT FOR YOUR PROJECT: having the two conflicting PC and CD
matrices had also caused ds9 to use the wrong WCS on the file. Open both the
'in.fits' and 'in-copy.fits' from the commands above in ds9 and look at the
WCS pointing directions of the file you attached and you'll see what I mean.

Finally, after you remove the PC matrix and do the conversion, you'll notice
that the PV coefficients will be different from what you initially posted from
the Python script! This is the Gnuastro TPV keywords:


PV1_1   =                  1.0 / TPV coefficient: x                           
 
PV1_4   =   6.719199510402E-02 / TPV coefficient: xx                          
 
PV1_5   =   4.672420249891E-02 / TPV coefficient: xy                          
 
PV1_6   =   4.679764407813E-02 / TPV coefficient: yy                          
 
PV2_1   =                  1.0 / TPV coefficient: y                           
 
PV2_4   =   3.174852324288E-02 / TPV coefficient: yy                          
 
PV2_5   =   7.274694257521E-02 / TPV coefficient: yx                          
 
PV2_6   =   3.601849423531E-02 / TPV coefficient: xx                    


And here are the ones you posted in the Python output:


PV1_1   =                  1.0                                                
 
PV1_4   =   0.1441040201614795                                                
 
PV1_5   =  0.02843416033032714                                                
 
PV1_6   =  0.03555357197501806                                                
 
PV2_1   =                  1.0                                                
 
PV2_4   =   0.0356690972763983                                                
 
PV2_5   =     0.07554122718085                                                
 
PV2_6   =  0.03865751488185146                                                
 


Since ds9 can correctly match the original image without the PC matrix and the
Gnuastro-TPV image, we can say that Gnuastro's conversion is correct. So
probably the Python script has also been confused by the conflicting PC and CD
matrices. 

You can get the most recent tarball that has this solution here. 

https://akhlaghi.org/gnuastro.pdf
https://akhlaghi.org/gnuastro-latest.tar.lz

About this new tarball, let me just say that NoiseChisel's thresholding and
sky subtraction step have been greatly improved in Commit 4a08b3fc205f
<http://git.savannah.gnu.org/cgit/gnuastro.git/commit/?id=4a08b3fc205f>. So if
you modified NoiseChisel's default parameters, make sure they are also useful
with this new version (it is much more robust, so I propose updating your
NoiseChisel scripts as soon as possible ;-)). I'll post a few examples in the
announcement for Gnuastro's 0.14 release.

    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?59459>

_______________________________________________
  Message sent via Savannah
  https://savannah.gnu.org/