[gnuastro-commits] master 524558e: Crop: --widthinpix option for width to be read as pixels in WCS mode

[Mohammad Akhlaghi]

branch: master
commit 524558e1a21d91144be86b8dc14fc9f0e91dd8b7
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    Crop: --widthinpix option for width to be read as pixels in WCS mode
    
    Until now, in Crop's WCS mode, the value to '--width' had to be in the same
    units as the WCS (degrees for RA, Dec). However, this was inconvenient when
    you had coordinates in RA/Dec, but wanted crops of a certain width (in
    pixels): you would have to get the pixel scale with 'astfits' for example,
    then convert your desired pixel-width to degrees.
    
    With this commit, Crop now has a new '--widthinpix' option (it is just an
    on-off option: meaning that it doesn't take a value). When this option is
    called, and Crop is in WCS mode, the value that the user gives to '--width'
    will be interpreted as pixels, not degrees.
    
    This new feature has also been included in the warning message that is
    printed when the crop size is larger than 50000 pixels in WCS mode. This
    usually happens when the user forgets that in WCS mode by default '--width'
    is in units of degrees, not pixels, and thus gives values like
    '--width=500'!
    
    This feature was suggested by Jesús Varela.
---
 NEWS                         |  6 +++
 THANKS                       |  1 +
 bin/crop/args.h              | 13 ++++++
 bin/crop/main.h              |  1 +
 bin/crop/ui.c                | 96 +++++++++++++++++++++++++++-----------------
 bin/crop/ui.h                |  3 +-
 doc/announce-acknowledge.txt |  1 +
 doc/gnuastro.texi            | 11 ++++-
 8 files changed, 93 insertions(+), 39 deletions(-)

diff --git a/NEWS b/NEWS
index a219776..3de3812 100644
--- a/NEWS
+++ b/NEWS
@@ -7,6 +7,12 @@ See the end of the file for license conditions.
 
 ** New features
 
+  Crop:
+   --widthinpix: when in WCS-mode, the value given to '--width' will be
+     interpreted in units of pixels. This will allow having crops of a
+     fixed pixel size, even when the mode of the central coordinate is in
+     WCS. This was suggested by Jesús Varela.
+
 ** Removed features
 
 ** Changed features
diff --git a/THANKS b/THANKS
index 6a28632..ef84352 100644
--- a/THANKS
+++ b/THANKS
@@ -109,6 +109,7 @@ support in Gnuastro. The list is ordered alphabetically (by 
family name).
     Ignacio Trujillo                     trujillo@iac.es
     Peter Teuben                         teuben@umd.edu
     David Valls-Gabaud                   david.valls-gabaud@obspm.fr
+    Jesús Varela                         jvarela@cefca.es
     Aaron Watkins                        aaron.watkins@oulu.fi
     Richard Wilbur                       richard.wilbur@gmail.com
     Michael H.F. Wilkinson               m.h.f.wilkinson@gmail.com
diff --git a/bin/crop/args.h b/bin/crop/args.h
index e820275..e4cf456 100644
--- a/bin/crop/args.h
+++ b/bin/crop/args.h
@@ -167,6 +167,19 @@ struct argp_option program_options[] =
       gal_options_parse_csv_float64
     },
     {
+      "widthinpix",
+      UI_KEY_WIDTHINPIX,
+      0,
+      0,
+      "--width is in pixels (even in WCS mode).",
+      GAL_OPTIONS_GROUP_OUTPUT,
+      &p->widthinpix,
+      GAL_OPTIONS_NO_ARG_TYPE,
+      GAL_OPTIONS_RANGE_0_OR_1,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
+    },
+    {
       "center",
       UI_KEY_CENTER,
       "FLT[,...]",
diff --git a/bin/crop/main.h b/bin/crop/main.h
index f18a88c..2e3016f 100644
--- a/bin/crop/main.h
+++ b/bin/crop/main.h
@@ -91,6 +91,7 @@ struct cropparams
   gal_data_t    *incheckcenter;  /* Value given to '--checkcenter'.       */
   gal_data_t           *center;  /* Center position of crop.              */
   gal_data_t            *width;  /* Width of crop when defined by center. */
+  uint8_t           widthinpix;  /* Given width is in pixels (in WCS mode)*/
   char                *catname;  /* Name of input catalog.                */
   char                 *cathdu;  /* HDU of catalog if its a FITS file.    */
   char                *namecol;  /* Filename (without suffix) of crop col.*/
diff --git a/bin/crop/ui.c b/bin/crop/ui.c
index 892e101..aefd16e 100644
--- a/bin/crop/ui.c
+++ b/bin/crop/ui.c
@@ -528,6 +528,24 @@ ui_check_options_and_arguments(struct cropparams *p)
 /**************************************************************/
 /***************       Preparations         *******************/
 /**************************************************************/
+#define UI_WIDTH_TOO_LARGE_SIZE 50000
+static void
+ui_warning_width_is_too_large(double width, size_t dim_num,
+                              double pixwidth, double pixscale)
+{
+  error(EXIT_SUCCESS, 0, "WARNING: value %g (requested WCS-based "
+        "width along dimension %zu) translates to %.0f pixels on "
+        "this dataset! This is most probably not what you wanted! "
+        "Note that the dataset's pixel size in this dimension is "
+        "%g. If you intended this number to show the width in "
+        "pixels, please add the '--widthinpix' option", width,
+        dim_num, pixwidth, pixscale);
+}
+
+
+
+
+
 /* When the crop is defined by its center, the final width that we need
    must be in actual number of pixels (an integer). But the user's values
    can be in WCS mode or even in image mode, they may be non-integers. */
@@ -574,32 +592,25 @@ ui_set_img_sizes(struct cropparams *p)
   /* WCS mode. */
   if(p->mode==IMGCROP_MODE_WCS)
     {
-      /* Fill in the widths depending on the mode. */
-      for(i=0;i<ndim;++i)
-        {
-          /* Convert the width in units of the input's WCS into pixels. */
-          pwidth = warray[i]/p->pixscale[i];
-          if(pwidth<1 || pwidth>50000)
-            error(EXIT_FAILURE, 0, "value %g (requested width along "
-                  "dimension %zu) translates to %.0f pixels on this "
-                  "dataset. This is probably not what you wanted. Note "
-                  "that the dataset's resolution in this dimension is "
-                  "%g.\n\n"
-                  "You can do the conversion to the dataset's WCS units "
-                  "prior to calling Crop. Alternatively, you can specify "
-                  "all the coordinates/sizes in image (not WCS) units and "
-                  "use the '--mode=img' option", warray[i], i+1, pwidth,
-                  p->pixscale[i]);
-
-          /* Write the single valued width in WCS and the image width for
-             this dimension. */
-          p->iwidth[i]=GAL_DIMENSION_FLT_TO_INT(pwidth);
-          if(p->iwidth[i]%2==0)
-            {
-              p->iwidth[i] += 1;
-              warray[i]    += p->pixscale[i];
-            }
-        }
+      /* Only convert the width if it is actually in degrees. */
+      if(p->widthinpix==0)
+        for(i=0;i<ndim;++i)
+          {
+            /* Calculate the pixel width. */
+            pwidth = warray[i]/p->pixscale[i];
+            if(pwidth<1 || pwidth>UI_WIDTH_TOO_LARGE_SIZE)
+              ui_warning_width_is_too_large(warray[i], i+1, pwidth,
+                                          p->pixscale[i]);
+
+            /* Write the single valued width in WCS and the image width for
+               this dimension. */
+            p->iwidth[i]=GAL_DIMENSION_FLT_TO_INT(pwidth);
+            if(p->iwidth[i]%2==0)
+              {
+                p->iwidth[i] += 1;
+                warray[i]    += p->pixscale[i];
+              }
+          }
 
       /* Checkcenter: */
       if(p->incheckcenter)
@@ -733,10 +744,10 @@ ui_read_cols(struct cropparams *p)
   /* Make sure more columns were not read (the name matchings might result
      in more than one column being read from the inputs). */
   if( gal_list_data_number(cols) != ndim + (p->namecol!=NULL) )
-    gal_tableintern_error_col_selection(p->catname, p->cathdu, "too many "
-                                        "columns were selected by the given "
-                                        "values to the options ending in "
-                                        "'col'.");
+    gal_tableintern_error_col_selection(p->catname, p->cathdu, "too "
+                                        "many columns were selected "
+                                        "by the given values to the "
+                                        "options ending in 'col'.");
 
 
   /* Put the information in each column in the proper place. */
@@ -879,7 +890,7 @@ ui_preparations_to_img_mode(struct cropparams *p)
 {
   size_t i;
   int nwcs;
-  double *darr, *pixscale;
+  double *darr, pixwidth, *pixscale;
   struct wcsprm *wcs=gal_wcs_read(p->inputs->v, p->cp.hdu,
                                   p->cp.wcslinearmatrix,
                                   p->hstartwcs, p->hendwcs, &nwcs);
@@ -911,15 +922,28 @@ ui_preparations_to_img_mode(struct cropparams *p)
               "defined by their center (with '--center' or '--catalog') a "
               "width is necessary (using the '--width' option)");
 
-      /* Check the requested width and convert it to pixels. */
-      darr=p->width->array;
+      /* Check if the dataset actually has WCS! */
       if(wcs->naxis<p->width->size)
         error(EXIT_FAILURE, 0, "%s (hdu %s): its WCS has %d dimensions "
               "but %zu dimensions given to '--width'", p->inputs->v,
               p->cp.hdu, wcs->naxis, p->width->size);
-      pixscale=gal_wcs_pixel_scale(wcs);
-      for(i=0;i<p->width->size;++i) darr[i] /= pixscale[i];
-      free(pixscale);
+
+      /* If the width is in WCS units, convert it to pixels. */
+      if(p->widthinpix==0)
+        {
+          darr=p->width->array;
+          pixscale=gal_wcs_pixel_scale(wcs);
+          for(i=0;i<p->width->size;++i)
+            {
+              pixwidth = darr[i] / pixscale[i];
+              if(pixwidth>UI_WIDTH_TOO_LARGE_SIZE)
+                ui_warning_width_is_too_large(darr[i], i+1, pixwidth,
+                                              pixscale[i]);
+              else
+                darr[i]=pixscale[i];
+            }
+          free(pixscale);
+        }
     }
 
   /* Switch the mode and return. */
diff --git a/bin/crop/ui.h b/bin/crop/ui.h
index 9a34fa2..35a6f58 100644
--- a/bin/crop/ui.h
+++ b/bin/crop/ui.h
@@ -46,7 +46,7 @@ enum program_args_groups
 /* Available letters for short options:
 
    a d e f g i j k m r t u v y
-   A B E G H J L Q R W X Y
+   A B E G H J L Q R W Y
 */
 enum option_keys_enum
 {
@@ -60,6 +60,7 @@ enum option_keys_enum
   UI_KEY_ZEROISNOTBLANK = 'z',
   UI_KEY_MODE           = 'O',
   UI_KEY_WIDTH          = 'w',
+  UI_KEY_WIDTHINPIX     = 'X',
   UI_KEY_CENTER         = 'c',
   UI_KEY_COORDCOL       = 'x',
 
diff --git a/doc/announce-acknowledge.txt b/doc/announce-acknowledge.txt
index edd0a87..f612c82 100644
--- a/doc/announce-acknowledge.txt
+++ b/doc/announce-acknowledge.txt
@@ -4,6 +4,7 @@ Sepideh Eskandarlou
 Zahra Hosseini
 Sebastian Luna-Valero
 Peter Teuben
+Jesús Varela
 
 
 
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 24010ad..1f7cbe2 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -12133,7 +12133,8 @@ Please see the description of this option in 
@ref{Invoking astcrop} for its synt
 @end table
 
 @item WCS coordinates
-In WCS mode (@option{--mode=wcs}), the coordinates and widths are interpreted 
using the World Coordinate System (WCS, that must accompany the dataset), not 
pixel coordinates.
+In WCS mode (@option{--mode=wcs}), the coordinates and width are interpreted 
using the World Coordinate System (WCS, that must accompany the dataset), not 
pixel coordinates.
+You can optionally use @option{--widthinpix} for the width to be interpreted 
in pixels (even though the coordinates are in WCS).
 In WCS mode, Crop accepts multiple datasets as input.
 When the cropped region (defined by its center or vertices) overlaps with 
multiple of the input images/tiles, the overlapping regions will be taken from 
the respective input (they will be stitched when necessary for each output 
crop).
 
@@ -12152,6 +12153,7 @@ There are multiple ways to define your cropped region 
as listed below.
 Similar to catalog inputs in Image mode (above), except that the values along 
each dimension are assumed to have the same units as the dataset's WCS 
information.
 For example, the central RA and Dec value for each crop will be read from the 
first and second calls to the @option{--coordcol} option.
 The width of the cropped box (in units of the WCS, or degrees in RA and Dec 
mode) must be specified with the @option{--width} option.
+You can optionally use @option{--widthinpix} for the value of @option{--width} 
to be interpreted in pixels.
 
 @item Center of a single crop (on the command-line)
 You can specify the center of only one crop box with the @option{--center} 
option.
@@ -12359,7 +12361,7 @@ The units of the coordinates are read based on the 
value to the @option{--mode}
 @item -w FLT[,FLT[,...]]
 @itemx --width=FLT[,FLT[,...]]
 Width of the cropped region about coordinate given to @option{--center}.
-If in WCS mode, value(s) given to this option will be read in the same units 
as the dataset's WCS information along this dimension.
+If in WCS mode, value(s) given to this option will be read in the same units 
as the dataset's WCS information along this dimension (unless 
@option{--widthinpix} is given).
 This option may take either a single value (to be used for all dimensions: 
@option{--width=10} in image-mode will crop a @mymath{10\times10} pixel image) 
or multiple values (a specific value for each dimension: @option{--width=10,20} 
in image-mode will crop a @mymath{10\times20} pixel image).
 
 The @code{--width} option also accepts fractions.
@@ -12372,6 +12374,11 @@ The basic reason for making an odd-sided crop is that 
your given central coordin
 When the crop has an odd number of pixels in each dimension, that pixel can be 
very well defined as the ``central'' pixel of the crop, making it unambiguously 
easy to identify.
 However, for an even-sided crop, it will be very hard to identify the central 
pixel (it can be on any of the four pixels adjacent to the central point of the 
image!).
 
+@item -X
+@itemx --widthinpix
+In WCS mode, interpret the value to @option{--width} as number of pixels, not 
the WCS units like degrees.
+This is useful when you want a fixed crop size in pixels, even though your 
center coordinates are in WCS (for example RA and Dec).
+
 @item -l STR
 @itemx -l FLT:FLT,...
 @itemx --polygon=STR