[bug #58974] WCS conversion not reasonable on processed TPV data

[Mohammad Akhlaghi]

URL:
  <https://savannah.gnu.org/bugs/?58974>

                 Summary: WCS conversion not reasonable on processed TPV data
                 Project: GNU Astronomy Utilities
            Submitted by: makhlaghi
            Submitted on: Tue 18 Aug 2020 02:36:32 AM BST
                Category: Libraries
                Severity: 3 - Normal
              Item Group: Output not reasonable
                  Status: In Progress
                 Privacy: Public
             Assigned to: makhlaghi
             Open/Closed: Open
         Discussion Lock: Any

    _______________________________________________________

Details:

When WCS headers with TPV distortions like below are used, Gnuastro won't
convert pixel to WCS coordinates (and vice-versa) properly! Producing
incorrect results. 


                      / World Coordinate System (WCS)                         
 
WCSAXES =                    2 / Number of coordinate axes                    
 
CRPIX1  =        1763.09429932 / Pixel coordinate of reference point          
 
CRPIX2  =         612.79893494 / Pixel coordinate of reference point          
 
PC1_1   =   -0.015946182574045 / Coordinate transformation matrix element     
 
PC1_2   =    0.999917332233843 / Coordinate transformation matrix element     
 
PC2_1   =   -0.999829574614546 / Coordinate transformation matrix element     
 
PC2_2   =   -0.015870824365013 / Coordinate transformation matrix element     
 
CDELT1  =  0.00078807879448870 / [deg] Coordinate increment at reference point
 
CDELT2  =  0.00078813807299668 / [deg] Coordinate increment at reference point
 
CUNIT1  = 'deg'                / Units of coordinate increment and value      
 
CUNIT2  = 'deg'                / Units of coordinate increment and value      
 
CTYPE1  = 'RA---TPV'           / TAN (gnomonic) projection + distortions      
 
CTYPE2  = 'DEC--TPV'           / TAN (gnomonic) projection + distortions      
 
CRVAL1  =       113.3064285324 / [deg] Coordinate value at reference point    
 
CRVAL2  =       31.87889767116 / [deg] Coordinate value at reference point    
 
LONPOLE =                180.0 / [deg] Native longitude of celestial pole     
 
LATPOLE =        31.8788976712 / [deg] Native latitude of celestial pole      
 
MJDREF  =                  0.0 / [d] MJD of fiducial time                     
 
RADESYS = 'ICRS'               / Equatorial coordinate system                 
 
PV1_0   =   9.714784609182E-05 / TPV coefficient: 1                           
 
PV1_1   =   9.999868977236E-01 / TPV coefficient: x                           
 
PV1_2   =  -8.842015406255E-06 / TPV coefficient: y                           
 
PV1_4   =   1.822783853090E-04 / TPV coefficient: xx                          
 
PV1_5   =   1.365776728868E-04 / TPV coefficient: xy                          
 
PV1_6   =   2.257022520459E-05 / TPV coefficient: yy                          
 
PV1_7   =   1.920222380359E-04 / TPV coefficient: xxx                         
 
PV1_8   =  -3.830056501464E-05 / TPV coefficient: xxy                         
 
PV1_9   =  -1.347022772508E-04 / TPV coefficient: xyy                         
 
PV1_10  =  -2.471225520868E-05 / TPV coefficient: yyy                         
 
PV2_0   =   4.672062523823E-04 / TPV coefficient: 1                           
 
PV2_1   =   1.000878428319E+00 / TPV coefficient: y                           
 
PV2_2   =  -4.973468945505E-04 / TPV coefficient: x                           
 
PV2_4   =  -1.411294363997E-04 / TPV coefficient: yy                          
 
PV2_5   =   7.439750263035E-04 / TPV coefficient: yx                          
 
PV2_6   =  -1.224454154415E-03 / TPV coefficient: xx                          
 
PV2_7   =  -7.014289760875E-04 / TPV coefficient: yyy                         
 
PV2_8   =   3.351859178423E-04 / TPV coefficient: yyx                         
 
PV2_9   =  -4.470842932692E-04 / TPV coefficient: yxx                         
 
PV2_10  =   9.945076014402E-04 / TPV coefficient: xxx 


After contacting the WCSLIB developer (Mark Calabretta), we found out that
this is actually a bug there: "the problem was that TPV assumes that CDi_ja is
used to define the linear transformation.  (Consequently, the independent
variables of the TPV polynomial are intermediate world coordinates rather than
intermediate pixel coordinates.  That's the wrong way to do it,but can be
excused as an historical anomaly.)  Thus, unlike other sequent distortions, if
PCi_ja is used instead of CDi_ja, then CDELTja must be applied before the TPV
polynomial. WCSLIB wasn't doing that.". Mark informed us that this will be
fixed in WCSLIB 7.4.

But it will take time (maybe a few years!) for the larger community to
transition to WCSLIB 7.4 and later. So to avoid confusing situations happening
for people in this transition period, in Gnuastro, when TPV keywords are
given, we should save a CD matrix, not PC+CDELT matrices. 

I am now working on this.




    _______________________________________________________

Reply to this item at:

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

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