[gnuastro-commits] master f4a643a: MakeCatalog: new --forcereadstd for read STD image even if not needed

[Mohammad Akhlaghi]

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

    MakeCatalog: new --forcereadstd for read STD image even if not needed
    
    Until now, MakeCatalog would only read the standard deviation (STD) image
    when one of the input columns needed it. However, one of MakeCatalog's
    output metadata (the surface brightenss limit) also needs the STD image and
    as a result nothing was printed. This could cause confusions.
    
    With this commit, there is a new `--forcereadstd' option that will ensure
    that the STD image is read even if no column needs it. Also, when its not
    present, instead of not printing anything at all, it now fills the same two
    lines with a notice of why there is no surface brightness limit.
    
    Also, while testing on the sample SDSS images, I noticed that the 10^4
    scale that we have set for a warning on WCS is a little too strict, so I
    loosened it to 10^5.
    
    This was added after a fruitful discussion with Raúl Infante-Sainz.
---
 NEWS                      | 13 ++++++-
 bin/mkcatalog/args.h      | 13 +++++++
 bin/mkcatalog/main.h      |  1 +
 bin/mkcatalog/mkcatalog.c | 94 +++++++++++++++++++++++++++--------------------
 bin/mkcatalog/ui.c        | 10 +++--
 bin/mkcatalog/ui.h        |  1 +
 doc/gnuastro.texi         | 10 ++++-
 lib/wcs.c                 |  6 +--
 8 files changed, 101 insertions(+), 47 deletions(-)

diff --git a/NEWS b/NEWS
index ad4c8e8..ee91c0a 100644
--- a/NEWS
+++ b/NEWS
@@ -36,7 +36,11 @@ See the end of the file for license conditions.
      standard.
 
   MakeCatalog:
-   --sigmaclip: defines the sigma-clipping parameters for those columns.
+   --sigmaclip: define sigma-clipping parameters for the new `--sigclip-*'
+     columns.
+   --forcereadstd: Read the standard deviation image even if not needed by
+     any columns. This is useful when you want the surface brightness
+     limit, but don't need any error-related columns.
    New output columns:
      --sigclip-number: Number of sigma-clipped pixels in object/clump.
      --sigclip-median: Sigma-clipped median of pixels in object/clump.
@@ -96,6 +100,13 @@ See the end of the file for license conditions.
      the start of the option name, makes it easier to find in the help list
      and also to understand generally.
 
+  MakeCatalog:
+   - Until now, if no standard deviation image was requested, MakeCatalog
+     wouldn't include any surface brightness limit metadata in its
+     output. Now, those two lines are filled, but with a notice on the
+     cause (that there is no standard deviation image), and suggesting
+     solutions.
+
   NoiseChisel:
    - Until now, when NoiseChisel didn't detect any pixels, it just printed
      a message and wouldn't not make any output file. This was very
diff --git a/bin/mkcatalog/args.h b/bin/mkcatalog/args.h
index 3ef319b..9c382ef 100644
--- a/bin/mkcatalog/args.h
+++ b/bin/mkcatalog/args.h
@@ -163,6 +163,19 @@ struct argp_option program_options[] =
       GAL_OPTIONS_NOT_SET
     },
     {
+      "forcereadstd",
+      UI_KEY_FORCEREADSTD,
+      0,
+      0,
+      "Read STD even if no columns need it.",
+      GAL_OPTIONS_GROUP_INPUT,
+      &p->forcereadstd,
+      GAL_OPTIONS_NO_ARG_TYPE,
+      GAL_OPTIONS_RANGE_0_OR_1,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
+    },
+    {
       "zeropoint",
       UI_KEY_ZEROPOINT,
       "FLT",
diff --git a/bin/mkcatalog/main.h b/bin/mkcatalog/main.h
index a85fda1..95ed39a 100644
--- a/bin/mkcatalog/main.h
+++ b/bin/mkcatalog/main.h
@@ -192,6 +192,7 @@ struct mkcatalogparams
   uint8_t         noclumpsort;  /* Don't sort the clumps catalog.       */
   float             zeropoint;  /* Zero-point magnitude of object.      */
   uint8_t            variance;  /* Input STD file is actually variance. */
+  uint8_t        forcereadstd;  /* Read STD even if not needed.         */
   uint8_t         subtractsky;  /* ==1: subtract the Sky from values.   */
   float           sfmagnsigma;  /* Surface brightness multiple of sigma.*/
   float             sfmagarea;  /* Surface brightness area (arcsec^2).  */
diff --git a/bin/mkcatalog/mkcatalog.c b/bin/mkcatalog/mkcatalog.c
index 55361ff..8c0bd86 100644
--- a/bin/mkcatalog/mkcatalog.c
+++ b/bin/mkcatalog/mkcatalog.c
@@ -462,6 +462,7 @@ mkcatalog_outputs_same_start(struct mkcatalogparams *p, int 
o0c1,
       gal_list_str_add(&comments, str, 0);
     }
 
+  /* Pixel area. */
   if(p->objects->wcs)
     {
       pixarea=gal_wcs_pixel_area_arcsec2(p->objects->wcs);
@@ -473,6 +474,7 @@ mkcatalog_outputs_same_start(struct mkcatalogparams *p, int 
o0c1,
         }
     }
 
+  /* Zeropoint magnitude */
   if(p->hasmag)
     {
       if( asprintf(&str, "Zeropoint magnitude: %.4f", p->zeropoint)<0 )
@@ -481,49 +483,53 @@ mkcatalog_outputs_same_start(struct mkcatalogparams *p, 
int o0c1,
     }
 
   /* Print surface brightness limits. */
-  if( !isnan(p->medstd) && !isnan(p->zeropoint) &&  !isnan(p->sfmagnsigma) )
+  if( !isnan(p->medstd) && !isnan(p->sfmagnsigma) )
     {
-      /* Per pixel. */
-      if( asprintf(&str, "%g sigma surface brightness (magnitude/pixel): "
-                   "%.3f", p->sfmagnsigma, ( -2.5f
-                                             *log10( p->sfmagnsigma
-                                                     * p->medstd )
-                                             + p->zeropoint ) )<0 )
-        error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
-      gal_list_str_add(&comments, str, 0);
-
-      /* Requested projected area: if a pixel area could be measured (a WCS
-         was given), then also estimate the surface brightness over one
-         arcsecond^2. From the pixel area, we know how many pixels are
-         necessary to fill the requested projected area (in
-         arcsecond^2). We also know that as the number of samples (pixels)
-         increases (to N), the noise increases by sqrt(N), see the full
-         discussion in the book. */
-      if(!isnan(pixarea) && !isnan(p->sfmagarea))
+      /* Only print magnitudes if a zeropoint is given. */
+      if( !isnan(p->zeropoint) )
         {
-          /* Prepare the comment/information. */
-          if(p->sfmagarea==1.0f)
-            tstr=NULL;
-          else
-            if( asprintf(&tstr, "%g-", p->sfmagarea)<0 )
-              error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
-          if( asprintf(&str, "%g sigma surface brightness "
-                       "(magnitude/%sarcsec^2): %.3f", p->sfmagnsigma,
-                       tstr ? tstr : "",
-                       ( -2.5f * log10( p->sfmagnsigma
-                                        * p->medstd
-                                        * sqrt( p->sfmagarea / pixarea) )
-                         + p->zeropoint ) )<0 )
+          /* Per pixel. */
+          if( asprintf(&str, "%g sigma surface brightness (magnitude/pixel): "
+                       "%.3f", p->sfmagnsigma, ( -2.5f
+                                                 *log10( p->sfmagnsigma
+                                                         * p->medstd )
+                                                 + p->zeropoint ) )<0 )
             error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
-
-          /* Add the final string/line to the catalog comments. */
           gal_list_str_add(&comments, str, 0);
 
-          /* Clean up (if necessary). */
-          if (tstr)
+          /* Requested projected area: if a pixel area could be measured (a
+             WCS was given), then also estimate the surface brightness over
+             one arcsecond^2. From the pixel area, we know how many pixels
+             are necessary to fill the requested projected area (in
+             arcsecond^2). We also know that as the number of samples
+             (pixels) increases (to N), the noise increases by sqrt(N), see
+             the full discussion in the book. */
+          if(!isnan(pixarea) && !isnan(p->sfmagarea))
             {
-              free(tstr);
-              tstr=NULL;
+              /* Prepare the comment/information. */
+              if(p->sfmagarea==1.0f)
+                tstr=NULL;
+              else
+                if( asprintf(&tstr, "%g-", p->sfmagarea)<0 )
+                  error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+              if( asprintf(&str, "%g sigma surface brightness "
+                           "(magnitude/%sarcsec^2): %.3f", p->sfmagnsigma,
+                           tstr ? tstr : "",
+                           ( -2.5f * log10( p->sfmagnsigma
+                                            * p->medstd
+                                            * sqrt( p->sfmagarea / pixarea) )
+                             + p->zeropoint ) )<0 )
+                error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+
+              /* Add the final string/line to the catalog comments. */
+              gal_list_str_add(&comments, str, 0);
+
+              /* Clean up (if necessary). */
+              if (tstr)
+                {
+                  free(tstr);
+                  tstr=NULL;
+                }
             }
         }
 
@@ -534,7 +540,17 @@ mkcatalog_outputs_same_start(struct mkcatalogparams *p, 
int o0c1,
         error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
       gal_list_str_add(&comments, str, 0);
     }
+  else
+    {
+      gal_checkset_allocate_copy("No surface brightness calcuations "
+                                 "because no STD image used.", &str);
+      gal_list_str_add(&comments, str, 0);
+      gal_checkset_allocate_copy("Ask for column that uses the STD image, "
+                                 "or `--forcereadstd'.", &str);
+      gal_list_str_add(&comments, str, 0);
+    }
 
+  /* The count-per-second correction. */
   if(p->cpscorr>1.0f)
     {
       if( asprintf(&str, "Counts-per-second correction: %.3f", p->cpscorr)<0 )
@@ -542,11 +558,11 @@ mkcatalog_outputs_same_start(struct mkcatalogparams *p, 
int o0c1,
       gal_list_str_add(&comments, str, 0);
     }
 
+  /* Print upper-limit parameters. */
   if(p->upperlimit)
     upperlimit_write_comments(p, &comments, 1);
 
-
-
+  /* Start column metadata. */
   if(p->cp.tableformat==GAL_TABLE_FORMAT_TXT)
     {
       if( asprintf(&str, "--------- Table columns ---------")<0 )
diff --git a/bin/mkcatalog/ui.c b/bin/mkcatalog/ui.c
index f8df810..01a54b7 100644
--- a/bin/mkcatalog/ui.c
+++ b/bin/mkcatalog/ui.c
@@ -913,8 +913,10 @@ ui_necessary_inputs(struct mkcatalogparams *p, int 
*values, int *sky,
 {
   size_t i;
 
+  /* Set necessary inputs based on options. */
+  if(p->forcereadstd) *std=1;
   if(p->upperlimit) *values=1;
-  if(p->spectrum)   *values=*std=1;
+  if(p->spectrum) *values=*std=1;
 
   /* Go over all the object columns. Note that the objects and clumps (if
      the `--clumpcat' option is given) inputs are mandatory and it is not
@@ -1177,10 +1179,12 @@ ui_preparations_read_inputs(struct mkcatalogparams *p)
         {
           for(column=p->objectcols; column!=NULL; column=column->next)
             if( !strcmp(column->unit, MKCATALOG_NO_UNIT) )
-              { free(column->unit); 
gal_checkset_allocate_copy(p->values->unit, &column->unit); }
+              { free(column->unit);
+                gal_checkset_allocate_copy(p->values->unit, &column->unit); }
           for(column=p->clumpcols; column!=NULL; column=column->next)
             if( !strcmp(column->unit, MKCATALOG_NO_UNIT) )
-              { free(column->unit); 
gal_checkset_allocate_copy(p->values->unit, &column->unit); }
+              { free(column->unit);
+                gal_checkset_allocate_copy(p->values->unit, &column->unit); }
         }
     }
 
diff --git a/bin/mkcatalog/ui.h b/bin/mkcatalog/ui.h
index bee10ba..0a6a33b 100644
--- a/bin/mkcatalog/ui.h
+++ b/bin/mkcatalog/ui.h
@@ -86,6 +86,7 @@ enum option_keys_enum
   UI_KEY_SKYHDU,
   UI_KEY_STDHDU,
   UI_KEY_WITHCLUMPS,
+  UI_KEY_FORCEREADSTD,
   UI_KEY_ZEROPOINT,
   UI_KEY_SIGMACLIP,
   UI_KEY_VARIANCE,
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index a45d009..2d1b8cd 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -14974,6 +14974,10 @@ The HDU of the Sky value standard deviation image.
 @item --variance
 The dataset given to @option{--stdfile} (and @option{--stdhdu} has the Sky 
variance of every pixel, not the Sky standard deviation.
 
+@item --forcereadstd
+Read the input STD image even if it is not required by any of the requested 
columns.
+This is because some of the output catalog's metadata may need it, for example 
to calculate the dataset's surface brightness limit (see @ref{Quantifying 
measurement limits}, configured with @option{--sfmagarea} and 
@option{--sfmagnsigma} in @ref{MakeCatalog output}).
+
 @item -z FLT
 @itemx --zeropoint=FLT
 The zero point magnitude for the input image, see @ref{Flux Brightness and 
magnitude}.
@@ -15525,11 +15529,15 @@ $ awk '!/^#/' out_c.txt | sort -g -k1,1 -k2,2
 @end example
 
 @item --sfmagnsigma=FLT
-The median standard deviation (from a @command{MEDSTD} keyword in the Sky 
standard deviation image) will be multiplied by the value to this option and 
its magnitude will be reported in the comments of the output catalog.
+Value to multiply with the median standard deviation (from a @command{MEDSTD} 
keyword in the Sky standard deviation image) for estimating the surface 
brightenss limit.
+Note that the surface brightness limit is only reported when a standard 
deviation image is read, in other words a column using it is requested (for 
example @option{--sn}) or @option{--forcereadstd} is called.
+
 This value is a per-pixel value, not per object/clump and is not found over an 
area or aperture, like the common @mymath{5\sigma} values that are commonly 
reported as a measure of depth or the upper-limit measurements (see 
@ref{Quantifying measurement limits}).
 
 @item --sfmagarea=FLT
 Area (in arcseconds squared) to convert the per-pixel estimation of 
@option{--sfmagnsigma} in the comments section of the output tables.
+Note that the surface brightness limit is only reported when a standard 
deviation image is read, in other words a column using it is requested (for 
example @option{--sn}) or @option{--forcereadstd} is called.
+
 Note that this is just a unit conversion using the World Coordinate System 
(WCS) information in the input's header.
 It does not actually do any measurements on this area.
 For random measurements on any area, please use the upper-limit columns of 
MakeCatalog (see the discussion on upper-limit measurements in @ref{Quantifying 
measurement limits}).
diff --git a/lib/wcs.c b/lib/wcs.c
index fb9676b..4eea2ce 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -747,12 +747,12 @@ gal_wcs_pixel_scale(struct wcsprm *wcs)
           }
 
       /* Do the check, print warning and make correction. */
-      if(maxrow!=minrow && maxrow/minrow>1e4 && warning_printed==0)
+      if(maxrow!=minrow && maxrow/minrow>1e5 && warning_printed==0)
         {
           fprintf(stderr, "\nWARNING: The input WCS matrix (possibly taken "
                   "from the FITS header keywords starting with `CD' or `PC') "
                   "contains values with very different scales (more than "
-                  "10^4 different). This is probably due to floating point "
+                  "10^5 different). This is probably due to floating point "
                   "errors. These values might bias the pixel scale (and "
                   "subsequent) calculations.\n\n"
                   "You can see the respective matrix with one of the "
@@ -766,7 +766,7 @@ gal_wcs_pixel_scale(struct wcsprm *wcs)
                   "You can delete the ones with obvious floating point "
                   "error values using the following command (assuming you "
                   "want to delete `CD1_2' and `CD2_1'). Afterwards, you can "
-                  "rerun your original command to remove this warning "
+                  "re-run your original command to remove this warning "
                   "message and possibly correct errors that it might have "
                   "caused.\n\n"
                   "    $ astfits file.fits --delete=CD1_2 --delete=CD2_1\n\n"