[gnuastro-commits] master 522f9a48: Library (wcs.h): wcsset is run after setting PC+CDELT

[Mohammad Akhlaghi]

branch: master
commit 522f9a489541b19951e74803a043f55a3fa1e6d6
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    Library (wcs.h): wcsset is run after setting PC+CDELT
    
    Until now, within 'gal_wcs_decompose_pc_cdelt' we would just update the
    'wcs->pc' and 'wcs->cdelt' values, but we weren't running 'wcsset'. As a
    result, WCSLIB was using some internal values prior to our change. This
    would result in in-correct outputs when the original rotation was defined
    by CROTA (which has been depreciated). For example in the following set of
    commands:
    
        # Download one (problematic!) dataset:
        raw=NGA_NGC1052-nd-int.fits.gz
        url1=http://galex.stsci.edu/data/GR6/pipe/01-vsn
        url2=05038-NGA_NGC1052/d/01-main/0001-img/07-try
        wget $url1/$url2/$img
    
        # Convert RA/Dec to X/Y (will be wrong!)
        echo 40.097200 -8.432200 \
             | asttable -c'arith $1 $2 wcs-to-img' \
                        --wcsfile=$img --wcshdu=0
    
        # Re-create the whole image (including WCS):
        astcrop $img -h0 --mode=img --section=:,: --output=good.fits
    
        # Convert RA/Dec to X/Y (will be correct!)
        echo 40.097200 -8.432200 \
             | asttable -c'arith $1 $2 wcs-to-img' \
                        --wcsfile=good.fits  --wcshdu=1
    
    After investigating the commands above, I noticed that the reason it gives
    a correct value on 'good.fits' is that we haven't been running 'wcssset'
    after correcting the internal structure information of the WCSLIB
    structure.
    
    With this commit, a call to 'wcsset' has been added after modifying the
    internal values in 'gal_wcs_decompose_pc_cdelt' and this fixed the problem.
    I also made two other minor edits:
    
     - In the book, there was no space added after the function name of
       'gal_wcs_write_in_fitsptr' (which is not the standard way we do it for
       all functions).
    
     - The error messages for 'wcsset' are now more standard (printing the name
       of the program, followed by function name, and then the WCSLIB message).
    
    This bug was reported by Ignacio Ruiz Cejudo and Raul Infante-Sainz.
    
    This fixes bug #62052.
---
 NEWS              | 10 ++++++----
 THANKS            |  1 +
 doc/gnuastro.texi |  2 +-
 lib/wcs.c         | 34 ++++++++++++++++++++--------------
 4 files changed, 28 insertions(+), 19 deletions(-)

diff --git a/NEWS b/NEWS
index 39122d2c..bb788f9f 100644
--- a/NEWS
+++ b/NEWS
@@ -186,15 +186,17 @@ See the end of the file for license conditions.
               in Julian dates for example) incorrectly read as float32
               (thus loosing precision); reported by Zohreh Ghaffari, fixed
               with help of Pedram Ashofteh Ardakani.
+  bug #61976: Table arithmetic date-to-sec produces same result for
+              different times (separated by about 1 second). This bug was
+              reported by Zohreh Ghaffari and Pedram Ashofteh Ardakani.
   bug #61967: DS9 polygon region files not read when they have width and
               color; reported by Zohreh Ghaffari, fixed by Pedram Ashofteh
               Ardakani.
   bug #62008: Arithmetic not deleting existing output when the 'makenew' is
               used (no FITS file exists in the reverse polish notation).
-  bug #61976: Table arithmetic date-to-sec produces same result for
-              different times (separated by about 1 second). This bug was
-              reported by Zohreh Ghaffari and Pedram Ashofteh Ardakani.
-
+  bug #62052: WCS decomposition of CD into PC+CDELT not setting interanal
+              values; reported by Ignacio Ruiz Cejudo and Raul
+              Infante-Sainz.
 
 
 
diff --git a/THANKS b/THANKS
index 4ecf01b6..72f7bd4e 100644
--- a/THANKS
+++ b/THANKS
@@ -94,6 +94,7 @@ support in Gnuastro. The list is ordered alphabetically (by 
family name).
     Bob Proulx                           bob@proulx.com
     Joseph Putko                         josephputko@gmail.com
     Samane Raji                          samaneraji@gmail.com
+    Ignacio Ruiz Cejudo                  igruiz04@ucm.es
     Francois Ochsenbein                  francois.ochsenbein@gmail.com
     Teymoor Saifollahi                   teymur.saif@gmail.com
     Joanna Sakowska                      js01093@surrey.ac.uk
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index fc1ccf4f..ada9b18e 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -29483,7 +29483,7 @@ When @code{extname!=NULL} it will be used as the FITS 
extension name.
 Any set of extra headers can also be written through the @code{headers} list 
and if @code{program_string!=NULL} it will be used in a commented keyword title 
just above the written version information.
 @end deftypefun
 
-@deftypefun void gal_wcs_write_in_fitsptr(fitsfile @code{*fptr}, struct wcsprm 
@code{*wcs})
+@deftypefun void gal_wcs_write_in_fitsptr (fitsfile @code{*fptr}, struct 
wcsprm @code{*wcs})
 Convert the input @code{wcs} structure (keeping the WCS programmatically) into 
FITS keywords and write them into the given FITS file pointer.
 This is a relatively low-level function which assumes the FITS file has 
already been opened with CFITSIO.
 If you just want to write the WCS into an empty file, you can use 
@code{gal_wcs_write} (which internally calls this function after creating the 
FITS file and later closes it safely).
diff --git a/lib/wcs.c b/lib/wcs.c
index cc7e0c64..4fd0c90d 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -220,11 +220,8 @@ gal_wcs_read_fitsptr(fitsfile *fptr, int linearmatrix, 
size_t hstartwcs,
   status=wcspih(fullheader, nkeys, relax, ctrl, &nreject, nwcs, &wcs);
   if(status)
     {
-      fprintf(stderr,"\n"
-              "##################\n"
-              "WCSLIB Warning: wcspih ERROR %d: %s.\n"
-              "##################\n",
-              status, wcs_errmsg[status]);
+      error(EXIT_SUCCESS, 0, "%s: WCSLIB Warning: wcspih ERROR %d: %s",
+            __func__, status, wcs_errmsg[status]);
       wcs=NULL; *nwcs=0;
     }
 
@@ -358,10 +355,8 @@ gal_wcs_read_fitsptr(fitsfile *fptr, int linearmatrix, 
size_t hstartwcs,
           status=wcsset(wcs);
           if(status)
             {
-              fprintf(stderr, "\n##################\n"
-                      "WCSLIB Warning: wcsset ERROR %d: %s.\n"
-                      "##################\n",
-                      status, wcs_errmsg[status]);
+              error(EXIT_SUCCESS, 0, "%s: WCSLIB warning: wcsset "
+                    "error %d: %s", __func__, status, wcs_errmsg[status]);
               wcsfree(wcs);
               wcs=NULL;
               *nwcs=0;
@@ -473,8 +468,8 @@ gal_wcs_create(double *crpix, double *crval, double *cdelt,
   /* Set up the wcs structure with the constants defined above. */
   status=wcsset(wcs);
   if(status)
-    error(EXIT_FAILURE, 0, "wcsset error %d: %s", status,
-          wcs_errmsg[status]);
+    error(EXIT_FAILURE, 0, "%s: wcsset error %d: %s", __func__,
+          status, wcs_errmsg[status]);
 
   /* If a CD matrix is desired make it. */
   if(linearmatrix==GAL_WCS_LINEAR_MATRIX_CD)
@@ -1620,6 +1615,11 @@ gal_wcs_warp_matrix(struct wcsprm *wcs)
       out[1] = -1 * wcs->cdelt[1] *sin(crota2);
       out[2] = wcs->cdelt[0] * sin(crota2);
       out[3] = wcs->cdelt[1] * cos(crota2);
+
+      /* For a check:
+      printf("cdelt: %f, %f\n", wcs->cdelt[0], wcs->cdelt[1]);
+      printf("cd:\n%f, %f\n%f, %f\n", out[0], out[1], out[2], out[3]);
+      */
     }
   else
     error(EXIT_FAILURE, 0, "%s: currently only PCi_ja and CDi_ja keywords "
@@ -1676,8 +1676,7 @@ gal_wcs_clean_errors(struct wcsprm *wcs)
 /* According to the FITS standard, in the 'PCi_j' WCS formalism, the matrix
    elements m_{ij} are encoded in the 'PCi_j' keywords and the scale
    factors are encoded in the 'CDELTi' keywords. There is also another
-   formalism (the 'CDi_j' formalism) which merges the two into one
-   matrix.
+   formalism (the 'CDi_j' formalism) which merges the two into one matrix.
 
    However, WCSLIB's internal operations are apparently done in the 'PCi_j'
    formalism. So its outputs are also all in that format by default. When
@@ -1689,6 +1688,7 @@ void
 gal_wcs_decompose_pc_cdelt(struct wcsprm *wcs)
 {
   size_t i, j;
+  int status=0;
   double *ps, *warp;
 
   /* If there is on WCS, then don't do anything. */
@@ -1700,7 +1700,7 @@ gal_wcs_decompose_pc_cdelt(struct wcsprm *wcs)
   if(ps==NULL) return;
 
   /* For a check.
-  printf("pc: %g, %g\n", pc[0], pc[1]);
+  printf("ps: %g, %g\n", ps[0], ps[1]);
   printf("warp: %g, %g, %g, %g\n", warp[0], warp[1],
          warp[2], warp[3]);
   */
@@ -1723,6 +1723,12 @@ gal_wcs_decompose_pc_cdelt(struct wcsprm *wcs)
      assumptioins about 'CDELTi'. */
   wcs->altlin=1;
 
+  /* Re-run 'wcsset' to update all the internal WCS parameters. */
+  status=wcsset(wcs);
+  if(status)
+    error(EXIT_SUCCESS, 0, "%s: WCSLIB warning: wcsset ERROR %d: %s",
+          __func__, status, wcs_errmsg[status]);
+
   /* Clean up. */
   free(ps);
   free(warp);