[gnuastro-commits] master ed41909 1/2: MakeCatalog: SBLMAG only added when pixel scale could be found

[Mohammad Akhlaghi]

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

    MakeCatalog: SBLMAG only added when pixel scale could be found
    
    Until now, if the pixel-scale couldn't be found (even when there was a WCS
    structure) MakeCatalog would crash (because the 'SBLMAG' keyword value
    would become NaN and CFITSIO doesn't write NaN values in numerical keywords
    (I guess this is a FITS convention). Such cases (where a WCS is present but
    the pixel scale can't be found) could happen in 3D datacubes for example
    because the 'gal_wcs_pixel_area_arcsec2' function would only return the
    pixel scale of 2D datasets.
    
    With this commit, the issue has been fixed in multiple layers one after
    another:
      1) To find the cause of the crash above, the 'gal_fits_key_write_in_ptr'
         function will now print a warning for any NaN valued keyword that is
         going to be written. A new small static function has been defined for
         this check.
      2) The 'gal_wcs_pixel_area_arcsec2' function will now accept datasets of
         any dimension, just as long as their first two dimensions are in units
         of degrees. The description of this function in the book is also
         corrected.
      3) In MakeCatalog, when we are writing the SBLMAG keyword, we explicitly
         check if a pixel-scale has been found or not (it may still happen for
         non RA/Dec images for example, like 2D histograms).
---
 bin/mkcatalog/mkcatalog.c | 19 +++++++++++--------
 doc/gnuastro.texi         |  8 ++++----
 lib/fits.c                | 34 ++++++++++++++++++++++++++++++++++
 lib/wcs.c                 |  7 +++----
 4 files changed, 52 insertions(+), 16 deletions(-)

diff --git a/bin/mkcatalog/mkcatalog.c b/bin/mkcatalog/mkcatalog.c
index 25d1f82..18f5a2d 100644
--- a/bin/mkcatalog/mkcatalog.c
+++ b/bin/mkcatalog/mkcatalog.c
@@ -501,14 +501,15 @@ mkcatalog_outputs_keys(struct mkcatalogparams *p, int 
o0c1)
                                          "Surface brightness limit per pixel.",
                                          "mag/pix");
 
-          /* Only print the SBL in fixed area if a WCS is present. */
-          if(p->objects->wcs)
+          /* Only print the SBL in fixed area if a WCS is present and a
+             pixel area could be deduced. */
+          if( !isnan(pixarea) )
             {
               /* Area used for measuring SBL. */
               mkcatalog_outputs_keys_numeric(&keylist, &p->sfmagarea,
                                              GAL_TYPE_FLOAT32, "SBLAREA",
-                                             "Area for surface brightness 
limit.",
-                                             "arcsec^2");
+                                             "Area for surface brightness "
+                                             "limit.", "arcsec^2");
 
               /* Per area, Surface brightness limit magnitude. */
               fvalue=gal_units_counts_to_mag(p->sfmagnsigma
@@ -523,10 +524,12 @@ mkcatalog_outputs_keys(struct mkcatalogparams *p, int 
o0c1)
             }
           else
             gal_fits_key_list_fullcomment_add_end(&keylist, "Can't "
-                   "write surface brightness limiting magnitude values "
-                   "in fixed area ('SBLAREA' and 'SBLMAG' keywords) "
-                   "because input doesn't have a world coordinate system "
-                   "to identify the pixel scale.", 0);
+                   "write surface brightness limiting magnitude (SBLM) "
+                   "values in fixed area ('SBLAREA' and 'SBLMAG' "
+                   "keywords) because input doesn't have a world "
+                   "coordinate system (WCS), or the first two "
+                   "coordinates of the WCS weren't angular positions "
+                   "in units of degrees.", 0);
         }
       else
         gal_fits_key_list_fullcomment_add_end(&keylist, "Can't write "
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 45ad0c6..4161e29 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -17065,7 +17065,7 @@ Plugging this into the magnitude equation, we get the 
@mymath{n\sigma} surface b
 @dispmath{SB_{{n\sigma,\rm A 
arcsec}^2}=-2.5\times\log_{10}{\left(n\sigma_m\over \sqrt{pA}\right)+z} 
\quad\quad [mag/arcsec^2]}
 
 @cindex World Coordinate System (WCS)
-MakeCatalog will calculate the input dataset's @mymath{SB_{n\sigma,\rm pixel}} 
and @mymath{SB_{{n\sigma,\rm A arcsec}^2}} and will write them as 
comments/meta-data in the output catalog(s)@footnote{Information on the general 
surface brightness limit is written in the @code{COMMENT} keywords if the 
requested output is a FITS table, or lines starting with a @code{#} when the 
output is plain-text.}, see @ref{Image surface brightness limit}.
+MakeCatalog will calculate the input dataset's @mymath{SB_{n\sigma,\rm pixel}} 
and @mymath{SB_{{n\sigma,\rm A arcsec}^2}} and will write them as the 
@code{SBLMAGPIX} and @code{SBLMAG} keywords the output catalog(s), see 
@ref{MakeCatalog output}.
 You can set your desired @mymath{n}-th multiple of @mymath{\sigma} and the 
@mymath{A} arcsec@mymath{^2} area using the following two options respectively: 
@option{--sfmagnsigma} and @option{--sfmagarea} (see @ref{MakeCatalog output}).
 Just note that @mymath{SB_{{n\sigma,\rm A arcsec}^2}} is only calculated if 
the input has World Coordinate System (WCS).
 Without WCS, the pixel scale cannot be derived.
@@ -25780,9 +25780,9 @@ return @code{NULL}.
 @end deftypefun
 
 @deftypefun double gal_wcs_pixel_area_arcsec2 (struct wcsprm @code{*wcs})
-Return the pixel area of @code{wcs} in arc-second squared. If the input WCS
-structure is not two dimensional and the units (@code{CUNIT} keywords) are
-not @code{deg} (for degrees), then this function will return a NaN.
+Return the pixel area of @code{wcs} in arc-second squared.
+This only works when the input dataset has atleast two dimensions and the 
units of the first two dimensions (@code{CUNIT} keywords) are @code{deg} (for 
degrees).
+In other cases, this function will return a NaN.
 @end deftypefun
 
 @deftypefun int gal_wcs_coverage (char @code{*filename}, char @code{*hdu}, 
size_t @code{*ondim}, double @code{**ocenter}, double @code{**owidth}, double 
@code{**omin}, double @code{**omax})
diff --git a/lib/fits.c b/lib/fits.c
index 646027a..0f31ec7 100644
--- a/lib/fits.c
+++ b/lib/fits.c
@@ -1899,6 +1899,36 @@ gal_fits_key_write(gal_fits_list_key_t **keylist, char 
*title,
 
 
 
+/* Fits doesn't allow NaN values, so if the type of float or double, we'll
+   just check to see if its NaN or not and let the user know the keyword
+   name (to help them fix it). */
+static void
+gal_fits_key_write_in_ptr_nan_check(gal_fits_list_key_t *tmp)
+{
+  int nanwarning=0;
+
+  /* Check the value. */
+  switch(tmp->type)
+    {
+    case GAL_TYPE_FLOAT32:
+      if( isnan( ((float *)(tmp->value))[0] ) ) nanwarning=1;
+      break;
+    case GAL_TYPE_FLOAT64:
+      if( isnan( ((double *)(tmp->value))[0] ) ) nanwarning=1;
+      break;
+    }
+
+  /* Print the warning. */
+  if(nanwarning)
+    error(EXIT_SUCCESS, 0, "%s: (WARNING) value of '%s' is NaN "
+          "and FITS doesn't recognize a NaN key value", __func__,
+          tmp->keyname);
+}
+
+
+
+
+
 /* Write the keywords in the gal_fits_list_key_t linked list to the FITS
    file. Every keyword that is written is freed, that is why we need the
    pointer to the linked list (to correct it after we finish). */
@@ -1928,6 +1958,10 @@ gal_fits_key_write_in_ptr(gal_fits_list_key_t **keylist, 
fitsfile *fptr)
           /* Write the basic key value and comments. */
           if(tmp->value)
             {
+              /* Print a warning if the value is NaN. */
+              gal_fits_key_write_in_ptr_nan_check(tmp);
+
+              /* Write/Update the keyword value. */
               if( fits_update_key(fptr, gal_fits_type_to_datatype(tmp->type),
                                   tmp->keyname, tmp->value, tmp->comment,
                                   &status) )
diff --git a/lib/wcs.c b/lib/wcs.c
index 0d2a508..f1d8144 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -1551,10 +1551,9 @@ gal_wcs_pixel_area_arcsec2(struct wcsprm *wcs)
   double out;
   double *pixscale;
 
-  /* A small sanity check. Later, when higher dimensions are necessary, we
-     can find which ones correlate to RA and Dec and use them to find the
-     pixel area in arcsec^2. */
-  if(wcs->naxis!=2) return NAN;
+  /* Some basic sanity checks. */
+  if(wcs==NULL) return NAN;
+  if(wcs->naxis==1) return NAN;
 
   /* Check if the units of the axis are degrees or not. Currently all FITS
      images I have worked with use 'deg' for degrees. If other alternatives