[gnuastro-commits] master c1eac22: Arithmetic: --wcsfile to specify file of output's WCS

[Mohammad Akhlaghi]

branch: master
commit c1eac228a9e4a12aa50c6b22d06977a335297a21
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>

    Arithmetic: --wcsfile to specify file of output's WCS
    
    With this new option, its now possible to specify a different file
    (possibly none of the inputs) to read the output's WCS from. It actually
    came up while trying to fix bug #55544 (described below).
    
    In the manual (version 0.8) we have mentioned that "output WCS information
    will be taken from the first input image encountered". This would logically
    mean the first dataset (with a WCS) that is found on the command-line.
    
    However, for implementation purposes (especially when there are hundreds of
    files), Arithmetic doesn't immediately read a file into memory! It will
    keep its reference and load it into memory only when processing is to be
    done on it (and immediately freed after it is no longer necessary).
    
    This caused confusion when the "where" operator is used, for example in a
    scenario like this:
    
    $ astarithmetic image1.fits --hdu=1 image2.fits --hdu=1 \
                    --output=masked.fits 1 eq nan where
    
    In this scenario, while the `image1.fits' is the first one on the
    command-line, `image2.fits' is the first one that is actually operated
    on. So Arithmetic would load the WCS of `image2.fits' and not `image1.fits'
    (as mentioned in the manual).
    
    To fix this bug, the reading of the output's WCS is now separated from the
    actual reading of the dataset arrays: it is now done when we add the files
    to the stack (without read the actual data into memory).
    
    This bug was reported by Raúl Infante-Sainz.
    
    This fixes bug #55544.
---
 NEWS                              |  6 ++++-
 bin/arithmetic/args.h             | 26 +++++++++++++++++++
 bin/arithmetic/arithmetic.c       |  2 +-
 bin/arithmetic/astarithmetic.conf |  3 +++
 bin/arithmetic/main.h             |  2 ++
 bin/arithmetic/operands.c         | 15 ++++++-----
 bin/arithmetic/ui.c               | 53 +++++++++++++++++++++++++++++++++++++++
 bin/arithmetic/ui.h               |  2 ++
 doc/gnuastro.texi                 | 46 ++++++++++++++++++++++++++-------
 9 files changed, 138 insertions(+), 17 deletions(-)

diff --git a/NEWS b/NEWS
index d2f3afe..b110e67 100644
--- a/NEWS
+++ b/NEWS
@@ -17,6 +17,10 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
      without changing the state of the operand stack (popping the top
      element). This can greatly help in debugging/understanding an
      Arithmetic command, especially as it gets longer, or more complicated.
+   --wcsfile and --wcshdu: these two options can be used to specify a
+     different file for reading the WCS of the output. This is useful when
+     the default (theh WCS of the first dataset that is read) is not the
+     required one.
 
   Fits:
    - Add "title" to group FITS keywords with `--write=/,"title name". This
@@ -65,7 +69,7 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
   bug #55439: Arithmetic segmentation fault when reusing dataset name.
   bug #55478: Memory mapping crashes when .gnuastro is not writable.
   bug #55491: NoiseChisel crash when no tiles good for quantile thresholding.
-
+  bug #55544: Arithmetic's output WCS with where operator is not as expected.
 
 
 
diff --git a/bin/arithmetic/args.h b/bin/arithmetic/args.h
index 1ef24dd..987148d 100644
--- a/bin/arithmetic/args.h
+++ b/bin/arithmetic/args.h
@@ -41,6 +41,32 @@ struct argp_option program_options[] =
       GAL_OPTIONS_NOT_MANDATORY,
       GAL_OPTIONS_NOT_SET
     },
+    {
+      "wcsfile",
+      UI_KEY_WCSFILE,
+      "STR",
+      0,
+      "File to use for output's WCS.",
+      GAL_OPTIONS_GROUP_INPUT,
+      &p->wcsfile,
+      GAL_TYPE_STRING,
+      GAL_OPTIONS_RANGE_ANY,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
+    },
+    {
+      "wcshdu",
+      UI_KEY_WCSHDU,
+      "STR",
+      0,
+      "HDU/extension to use for output's WCS.",
+      GAL_OPTIONS_GROUP_INPUT,
+      &p->wcshdu,
+      GAL_TYPE_STRING,
+      GAL_OPTIONS_RANGE_ANY,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
+    },
 
 
 
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index d3b5dc9..886d61c 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -1181,9 +1181,9 @@ reversepolish(struct arithmeticparams *p)
       filename=p->operands->filename;
       if( gal_fits_name_is_fits(filename) )
         {
+          /* Read the data, note that the WCS has already been set. */
           p->operands->data=gal_array_read_one_ch(filename, hdu, NULL,
                                                   p->cp.minmapsize);
-          p->refdata.wcs=gal_wcs_read(filename, hdu, 0, 0, &p->refdata.nwcs);
           if(!p->cp.quiet) printf(" - %s (hdu %s) is read.\n", filename, hdu);
         }
       else
diff --git a/bin/arithmetic/astarithmetic.conf 
b/bin/arithmetic/astarithmetic.conf
index 7b82ea6..ff83a14 100644
--- a/bin/arithmetic/astarithmetic.conf
+++ b/bin/arithmetic/astarithmetic.conf
@@ -16,3 +16,6 @@
 # permitted in any medium without royalty provided the copyright notice and
 # this notice are preserved.  This file is offered as-is, without any
 # warranty.
+
+# Inputs
+ wcshdu        1
\ No newline at end of file
diff --git a/bin/arithmetic/main.h b/bin/arithmetic/main.h
index 73110d8..d6f6661 100644
--- a/bin/arithmetic/main.h
+++ b/bin/arithmetic/main.h
@@ -73,6 +73,8 @@ struct arithmeticparams
   /* Input: */
   gal_list_str_t     *hdus;  /* List of all given HDU strings.          */
   gal_list_str_t   *tokens;  /* List of all arithmetic tokens.          */
+  char            *wcsfile;  /* File to use for output's WCS.           */
+  char             *wcshdu;  /* Extension to use for output's WCS.      */
   size_t        popcounter;  /* The number of FITS images popped.       */
   gal_data_t       refdata;  /* Container for information of the data.  */
   char          *globalhdu;  /* Single HDU for all inputs.              */
diff --git a/bin/arithmetic/operands.c b/bin/arithmetic/operands.c
index dc4f31c..a5071f7 100644
--- a/bin/arithmetic/operands.c
+++ b/bin/arithmetic/operands.c
@@ -304,6 +304,15 @@ operands_add(struct arithmeticparams *p, char *filename, 
gal_data_t *data)
                 gal_checkset_allocate_copy(p->globalhdu, &newnode->hdu);
               else
                 newnode->hdu=gal_list_str_pop(&p->hdus);
+
+              /* If no WCS is set yet, use the WCS of this image. */
+              if(p->refdata.wcs==NULL)
+                {
+                  p->refdata.wcs=gal_wcs_read(filename, newnode->hdu, 0, 0,
+                                              &p->refdata.nwcs);
+                  if(p->refdata.wcs && !p->cp.quiet)
+                    printf(" - WCS: %s (hdu %s).\n", filename, newnode->hdu);
+                }
             }
           else newnode->hdu=NULL;
         }
@@ -349,12 +358,6 @@ operands_pop(struct arithmeticparams *p, char *operator)
          identify variables. */
       if(data->name) { free(data->name); data->name=NULL; }
 
-      /* In case this is the first image that is read, then keep the WCS
-         information in the `refdata' structure.  */
-      if(p->popcounter==0)
-        p->refdata.wcs=gal_wcs_read(filename, hdu, 0, 0,
-                                    &p->refdata.nwcs);
-
       /* When the reference data structure's dimensionality is non-zero, it
          means that this is not the first image read. So, write its basic
          information into the reference data structure for future
diff --git a/bin/arithmetic/ui.c b/bin/arithmetic/ui.c
index 8b84588..c4baea8 100644
--- a/bin/arithmetic/ui.c
+++ b/bin/arithmetic/ui.c
@@ -28,6 +28,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <stdio.h>
 #include <string.h>
 
+#include <gnuastro/wcs.h>
 #include <gnuastro/list.h>
 #include <gnuastro/fits.h>
 #include <gnuastro/table.h>
@@ -225,6 +226,26 @@ parse_opt(int key, char *arg, struct argp_state *state)
 /**************************************************************/
 /***************       Sanity Check         *******************/
 /**************************************************************/
+static void
+ui_read_check_only_options(struct arithmeticparams *p)
+{
+  if(p->wcsfile)
+    {
+      if(gal_fits_name_is_fits(p->wcsfile)==0)
+        error(EXIT_FAILURE, 0, "%s: file given to `--wcsfile' must be in "
+              "FITS format with a recognizable FITS format suffix.",
+              p->wcsfile);
+      if(p->wcshdu==NULL)
+        error(EXIT_FAILURE, 0, "%s: no HDU/extension specified (file given "
+              "to `--wcsfile')! Please use `--wcshdu' to specify a "
+              "HDU/extension to read from", p->wcsfile);
+    }
+}
+
+
+
+
+
 /* Sanity check on options AND arguments. If only option values are to be
    checked, use `ui_read_check_only_options'. */
 static void
@@ -327,6 +348,30 @@ ui_check_options_and_arguments(struct arithmeticparams *p)
 
 
 
+static void
+ui_preparations(struct arithmeticparams *p)
+{
+  /* In case a file is specified to read the WCS from (and ignore input
+     datasets), read the WCS prior to starting parsing of the arguments. */
+  if(p->wcsfile)
+    {
+      p->refdata.wcs=gal_wcs_read(p->wcsfile, p->wcshdu, 0, 0,
+                                  &p->refdata.nwcs);
+      if(p->refdata.wcs)
+        {
+          if(!p->cp.quiet)
+            printf(" - WCS: %s (hdu %s).\n", p->wcsfile, p->wcshdu);
+        }
+      else
+        fprintf(stderr, "WARNING: %s (hdu %s) didn't contain a "
+                "(readable by WCSLIB) WCS.", p->wcsfile, p->wcshdu);
+    }
+}
+
+
+
+
+
 
 
 
@@ -392,10 +437,18 @@ ui_read_check_inputs_setup(int argc, char *argv[], struct 
arithmeticparams *p)
   gal_options_print_state(cp);
 
 
+  /* Sanity check only on options. */
+  ui_read_check_only_options(p);
+
+
   /* Check that the options and arguments fit well with each other. Note
      that arguments don't go in a configuration file. So this test should
      be done after (possibly) printing the option values. */
   ui_check_options_and_arguments(p);
+
+
+  /* Initial preparations. */
+  ui_preparations(p);
 }
 
 
diff --git a/bin/arithmetic/ui.h b/bin/arithmetic/ui.h
index 6ee011a..7263411 100644
--- a/bin/arithmetic/ui.h
+++ b/bin/arithmetic/ui.h
@@ -40,6 +40,8 @@ enum option_keys_enum
   /* With short-option version. */
   UI_KEY_GLOBALHDU       = 'g',
   UI_KEY_ONEDASIMAGE     = 'O',
+  UI_KEY_WCSFILE         = 'w',
+  UI_KEY_WCSHDU          = 'W'
 
   /* Only with long version (start with a value 1000, the rest will be set
      automatically). */
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 371ef36..7dc3584 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -12821,15 +12821,21 @@ $ astarithmetic img1.fits img2.fits img3.fits median  
              \
                 -h0 -h1 -h2 --out=median.fits
 @end example
 
-If the output is not a single number, and the @option{--output} option is
-not given, automatic output will use the name of the first FITS image
-encountered to generate an output file name, see @ref{Automatic
-output}. Also, output WCS information will be taken from the first input
-image encountered. When the output is a single number, that number will be
-printed in the standard output and no output file will be
-created. Arithmetic's notation for giving operands to operators is
-described in @ref{Reverse polish notation}. To ignore certain pixels, set
-them as blank, see @ref{Blank pixels}, for example with the @command{where}
+The output is last remaining operand on the stack, see @ref{Reverse polish
+notation}. When its a single number, it will be printed on the
+command-line. When the output is an array, it will be stored as a file. The
+name of the file can be specified with the @option{--output} option, but if
+its not given, Arithmetic will use ``automatic output'' on the name of the
+first FITS image encountered to generate an output file name, see
address@hidden output}.
+
+By default, the world coordinate system (WCS) information will be taken
+from the first input image (that contains a WCS) on the command-line. This
+can be modified with the @option{--wcsfile} option described below.
+
+Arithmetic's notation for giving operands to operators is fully described
+in @ref{Reverse polish notation}. To ignore certain pixels, set them as
+blank, see @ref{Blank pixels}, for example with the @command{where}
 operator (see @ref{Arithmetic operators}). See @ref{Common options} for a
 review of the options in all Gnuastro programs. Arithmetic just redefines
 the @option{--hdu} option as explained below.
@@ -12878,6 +12884,28 @@ interest is in the same HDU of all the files. When 
this option is called,
 any values given to the @option{--hdu} option (explained above) are ignored
 and will not be used.
 
address@hidden --wcsfile STR
+FITS Filename containing the WCS structure that must be written to the
+output. The HDU/extension should be specified with @option{--wcshdu}.
+
+When this option is used, the respective WCS will be read before any
+processing is done on the command-line and directly used in the final
+output. If the given file doesn't have any WCS, then the default WCS (first
+file on the command-line) will be used in the output.
+
+This option will mostly be used when the default file (first of the set of
+inputs) is not the one containing your desired WCS. But with this option,
+you can use Arithmetic to rewrite/chage the WCS of an existing FITS file
+from another file a command like this:
+
address@hidden
+$ astarithmetic data.fits --wcsfile=other.fits -ofinal.fits
address@hidden example
+
address@hidden --wcshdu STR
+HDU/extension to read the WCS within the file given to
address@hidden For more, see the description of @option{--wcsfile}.
+
 @item -O
 @itemx --onedasimage
 When final dataset to write as output only has one dimension, write it as a