[task #15676] Calculating SIP to TPV coversions

[Raul Infante-Sainz]

Follow-up Comment #8, task #15676 (project gnuastro):

I have checked the feature with a random image (only one image!) and there are
some things I would like to point out and let you know.

It works without crashing when reading the header from an image (header +
data). However, sometimes we only have the header with the WCS information as
input (for example if you obtain the WCS using `solve-field' program of
Astrometry.net). In that case, when reading from headers (without data,
attached file `header-sip.txt`), it crashed with this output:


$ astfits header-sip.fits -h0 --wcsdistortion=TPV

astfits: header-sip.fits (hdu: 0) has 0 dimensions! The most common cause for
this is a wrongly specified HDU. In some FITS images, the first HDU doesn't
have any data, the data is in subsequent extensions. So probably reading the
second HDU (with '--hdu=1' or '-h1') will solve the problem (following
CFITSIO's convention, currently HDU counting starts from 0).


It would be worth if it could read the data from headers only.


Now, let's use an image with data (header + data). Here it works without
crashing :-) The header of that image is attached in file `original-sip.txt'.


                      / World Coordinate System (WCS)
WCSAXES =                    2 / Number of coordinate axes
CRPIX1  =        976.816238403 / Pixel coordinate of reference point
CRPIX2  =        440.521743774 / Pixel coordinate of reference point
PC1_1   =    -1.00022174065751 / Coordinate transformation matrix element
PC1_2   =    -0.00319660295795 / Coordinate transformation matrix element
PC2_1   =    -0.00235639222511 / Coordinate transformation matrix element
PC2_2   =     0.99977077772003 / Coordinate transformation matrix element
CDELT1  =  9.34302943928898E-05 / [deg] Coordinate increment at reference
point
CDELT2  =  9.35240534339602E-05 / [deg] Coordinate increment at reference
point
CUNIT1  = 'deg'                / Units of coordinate increment and value
CUNIT2  = 'deg'                / Units of coordinate increment and value
CTYPE1  = 'RA---TAN-SIP'       / TAN (gnomonic) projection + SIP distortions
CTYPE2  = 'DEC--TAN-SIP'       / TAN (gnomonic) projection + SIP distortions
CRVAL1  =        39.9265010846 / [deg] Coordinate value at reference point
CRVAL2  =       -8.12618217247 / [deg] Coordinate value at reference point
LONPOLE =                180.0 / [deg] Native longitude of celestial pole
LATPOLE =       -8.12618217247 / [deg] Native latitude of celestial pole
RADESYS = 'FK5'                / Equatorial coordinate system
EQUINOX =               2000.0 / [yr] Equinox of equatorial coordinates
A_ORDER =                    2 / SIP polynomial degree, axis 1, pixel-to-sky
A_0_2   =    4.19888064518E-07 / SIP distortion coefficient: yy
A_1_1   =   -2.33927038690E-07 / SIP distortion coefficient: xy
A_2_0   =   -2.54499914902E-07 / SIP distortion coefficient: xx
AP_ORDER=                    2 / SIP polynomial degree, axis 1, sky-to-pixel
AP_0_0  =    2.98764311704E-05 / SIP inverse coefficient:
AP_0_1  =    1.67505488187E-07 / SIP inverse coefficient: y
AP_0_2  =   -4.20424622199E-07 / SIP inverse coefficient: yy
AP_1_0  =    9.45864856941E-08 / SIP inverse coefficient: x
AP_1_1  =    2.32274524960E-07 / SIP inverse coefficient: xy
AP_2_0  =    2.54581796288E-07 / SIP inverse coefficient: xx
B_ORDER =                    2 / SIP polynomial degree, axis 2, pixel-to-sky
B_0_2   =   -2.40565772303E-07 / SIP distortion coefficient: yy
B_1_1   =   -1.63905341261E-06 / SIP distortion coefficient: xy
B_2_0   =   -5.30068612640E-07 / SIP distortion coefficient: xx
BP_ORDER=                    2 / SIP polynomial degree, axis 2, sky-to-pixel
BP_0_0  =   -6.32530864314E-04 / SIP inverse coefficient:
BP_0_1  =    1.43043299047E-06 / SIP inverse coefficient: y
BP_0_2  =    2.39626050245E-07 / SIP inverse coefficient: yy
BP_1_0  =    6.82806205992E-07 / SIP inverse coefficient: x
BP_1_1  =    1.64047618237E-06 / SIP inverse coefficient: xy
BP_2_0  =    5.32170294564E-07 / SIP inverse coefficient: xx



This header comes from the following command:


astarithmetic image-no-wcs.fits --wcsfile=header-sip.fits --wcshdu=0
--output=original-sip.fits


With the previous command, I want to inject the WCS information of
`header-sip.fits' header obtained with Astrometry.net `solve-field' program
(so, `header-sip.txt' of attached files) into the image without the WCS
information. The result header of `original-sip.fits' image is what I attach
as `original-sip.txt'.

In short, I have the previous WCS/SIP parameters and now I transform them to
TPV.

With Python sip_tpv:


                      / World Coordinate System (WCS)
WCSAXES =                    2 / Number of coordinate axes
CRPIX1  =        976.816238403 / Pixel coordinate of reference point
CRPIX2  =        440.521743774 / Pixel coordinate of reference point
PC1_1   =    -1.00022174065751 / Coordinate transformation matrix element
PC1_2   =    -0.00319660295795 / Coordinate transformation matrix element
PC2_1   =    -0.00235639222511 / Coordinate transformation matrix element
PC2_2   =     0.99977077772003 / Coordinate transformation matrix element
CDELT1  =  9.34302943928898E-05 / [deg] Coordinate increment at reference
point
CDELT2  =  9.35240534339602E-05 / [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  =        39.9265010846 / [deg] Coordinate value at reference point
CRVAL2  =       -8.12618217247 / [deg] Coordinate value at reference point
LONPOLE =                180.0 / [deg] Native longitude of celestial pole
LATPOLE =       -8.12618217247 / [deg] Native latitude of celestial pole
RADESYS = 'FK5'                / Equatorial coordinate system
EQUINOX =               2000.0 / [yr] Equinox of equatorial coordinates
PV1_1   =                   1.0 / TPV coefficient: x
PV1_4   =  2.74744439572126E-03 / TPV coefficient: xx
PV1_5   = -2.51914579463917E-03 / TPV coefficient: xy
PV1_6   = -4.48803400336329E-03 / TPV coefficient: yy
PV2_1   =                   1.0 / TPV coefficient: y
PV2_4   = -2.52742435899630E-03 / TPV coefficient: yy
PV2_5   =  1.75088471295910E-02 / TPV coefficient: yx
PV2_6   = -5.71014067284307E-03 / TPV coefficient: xx


With Gnuastro --wcsdistortion=TPV:


                      / World Coordinate System (WCS)
WCSAXES =                    2 / Number of coordinate axes
CRPIX1  =         976.81623840 / Pixel coordinate of reference point
CRPIX2  =         440.52174377 / Pixel coordinate of reference point
PC1_1   =    -1.01912505239904 / Coordinate transformation matrix element
PC1_2   =    -0.75125324648638 / Coordinate transformation matrix element
PC2_1   =     0.02503243359873 / Coordinate transformation matrix element
PC2_2   =     0.99968663768025 / Coordinate transformation matrix element
CDELT1  =  2.93054865755093E-07 / [deg] Coordinate increment at reference
point
CDELT2  =  3.73525871020151E-03 / [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  =        39.9265010846 / [deg] Coordinate value at reference point
CRVAL2  =        -8.1261821725 / [deg] Coordinate value at reference point
LONPOLE =                180.0 / [deg] Native longitude of celestial pole
LATPOLE =        -8.1261821725 / [deg] Native latitude of celestial pole
RADESYS = 'FK5'                / Equatorial coordinate system
EQUINOX =               2000.0 / [yr] Equinox of equatorial coordinates
PV1_0   =  1.39396742233200E-10 / TPV coefficient: 1
PV1_1   =  1.00000050383900E+00 / TPV coefficient: x
PV1_2   = 7.61504572502400E-168 / TPV coefficient: y
PV1_4   = -4.62686729211900E-01 / TPV coefficient: xx
PV1_5   = -2.88130137773400E-161 / TPV coefficient: xy
PV1_6   = 1.97626258336499E-323 / TPV coefficient: yy
PV2_0   = 1.69617046769700E+151 / TPV coefficient: 1
PV2_1   =  1.00000092659400E+00 / TPV coefficient: y
PV2_2   = 6.13068244241300E+154 / TPV coefficient: x
PV2_4   = 2.53998312986400E-161 / TPV coefficient: yy
PV2_5   = -2.86314203367400E+00 / TPV coefficient: yx
PV2_6   = -1.59987324584300E+161 / TPV coefficient: xx


As you can see, there are some parameters which are very big/small!! And they
are not present in the PV parameters obtained with the Python script (PV1_2,
for example). What do you think? Let me know if you need more
information/inputs. Thanks a lot for the great work ;-)

(file #49343, file #49344, file #49345, file #49346)
    _______________________________________________________

Additional Item Attachment:

File name: header-sip.txt                 Size:5 KB
    <https://savannah.gnu.org/file/header-sip.txt?file_id=49343>

File name: python-pv.txt                  Size:1 KB
    <https://savannah.gnu.org/file/python-pv.txt?file_id=49344>

File name: original-sip.txt               Size:2 KB
    <https://savannah.gnu.org/file/original-sip.txt?file_id=49345>

File name: gnuastro-pv.txt                Size:2 KB
    <https://savannah.gnu.org/file/gnuastro-pv.txt?file_id=49346>



    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/task/?15676>

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