[gnuastro-commits] master 6653ac9 1/2: library (fits.h): new function to check if data is HEALPix

[Mohammad Akhlaghi]

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

    library (fits.h): new function to check if data is HEALPix
    
    2D HEALPix datasets (which have spherical "pixel"s and are used for
    whole-sky datasets like those of Plank) are stored as 1D arrays (mostly a
    table column). Since there isn't any internal solution in Gnuastro yet to
    convert 1D HEALPix datasets to a 2D image, their special storage would give
    un-expected errors for the users that could be confusing. Also, WCSLIB
    builds and installs a "utility" called 'HPXcvt' for this job, but most
    people are unaware of it.
    
    With this commit, a new function has been added to the Gnuastro library
    called 'gal_fits_hdu_is_healpix'. This function is used in the core library
    function to read FITS images and if the user's input is actually a 1D
    HEALPix dataset, it will print an informative error and abort. In the error
    message, it recomments WCSLIB's 'HPXcvt' program as a solution.
    
    Also, in the FITS program, when no action is requested (and it prints the
    basic information about every HDU), an extra comment column has been added
    that will check for HEALPix. If the dataset has the standard HEALPix
    keywords, the new comment column will mention this for the user to know.
    
    This issue was found with the help of Ignacio Trujillo.
---
 NEWS                      |  8 ++++++++
 bin/arithmetic/operands.c | 19 ++++++++++++++++---
 bin/fits/fits.c           | 34 ++++++++++++++++++++++++++++++----
 doc/gnuastro.texi         | 10 +++++++++-
 lib/fits.c                | 41 +++++++++++++++++++++++++++++++++++++++--
 lib/gnuastro/fits.h       |  3 +++
 6 files changed, 105 insertions(+), 10 deletions(-)

diff --git a/NEWS b/NEWS
index 0930039..d0a5d88 100644
--- a/NEWS
+++ b/NEWS
@@ -38,6 +38,9 @@ See the end of the file for license conditions.
      arithmetic:
          echo "113.64812416667 31.88732" \
               | asttable -c'arith $1 degree-to-ra $2 degree-to-dec'
+   - If input is a HEALpix grid (1D table column that represents the 2D
+     spherical representation of datasets), the programs will print a
+     warning, suggesting to use the 'HPXcvt' utility of WCSLIB.
 
   Arithmetic:
    - New operators:
@@ -113,6 +116,7 @@ See the end of the file for license conditions.
    - gal_blank_flag_remove: Remove all flagged elements in a dataset.
    - gal_blank_remove_rows: remove all rows that have at least one blank.
    - gal_dimension_dist_elliptical: Elliptical dist. of a point from center.
+   - gal_fits_hdu_is_healpix: Return 1 if HDU is a HEALpix grid.
    - gal_fits_hdu_datasum: calculate DATASUM of given HDU in given FITS file.
    - gal_fits_hdu_datasum_ptr: calculate DATASUM of opened FITS file pointer.
    - gal_pointer_allocate_ram_or_mmap: allocate space either in RAM or as a
@@ -143,6 +147,10 @@ See the end of the file for license conditions.
    - The '--pixelscale' option also prints the pixel area (for 2D inputs,
      or 2D slices of 3D inputs) and the voxel volume (for 3D inputs). Until
      now, it would only print the pixel scale along each dimension.
+   - When printing FITS file HDU information (no options given), a new
+     "Comments" column may be printed for each HDU in the end of the line.
+     It will be printed if special situations are found (for example a 2D
+     HEALPix grid, that is usually stored as a 1D array/column).
 
   NoiseChisel & Statistics:
    - New algorithm used to reject outlying tiles. In NoiseChisel this is
diff --git a/bin/arithmetic/operands.c b/bin/arithmetic/operands.c
index 9907c38..e019bc9 100644
--- a/bin/arithmetic/operands.c
+++ b/bin/arithmetic/operands.c
@@ -312,13 +312,26 @@ operands_add(struct arithmeticparams *p, char *filename, 
gal_data_t *data)
               readwcs = (p->wcsfile && !strcmp(p->wcsfile,"none")) ? 0 : 1;
               if(readwcs && p->refdata.wcs==NULL)
                 {
-                  dsize=gal_fits_img_info_dim(filename, newnode->hdu, &ndim);
+                  /* If the HDU is an image, read its size. */
+                  dsize = ( gal_fits_hdu_format(filename, 
newnode->hdu)==IMAGE_HDU
+                            ? gal_fits_img_info_dim(filename, newnode->hdu, 
&ndim)
+                            : NULL);
+
+                  /* Read the WCS. */
                   p->refdata.wcs=gal_wcs_read(filename, newnode->hdu, 0, 0,
                                               &p->refdata.nwcs);
-                  ndim=gal_dimension_remove_extra(ndim, dsize, p->refdata.wcs);
+
+                  /* Remove extra (length of 1) dimensions (if we had an
+                     image HDU). */
+                  if(dsize)
+                    {
+                      ndim=gal_dimension_remove_extra(ndim, dsize, 
p->refdata.wcs);
+                      free(dsize);
+                    }
+
+                  /* Let the user know that the WCS is read. */
                   if(p->refdata.wcs && !p->cp.quiet)
                     printf(" - WCS: %s (hdu %s).\n", filename, newnode->hdu);
-                  free(dsize);
                 }
             }
           else newnode->hdu=NULL;
diff --git a/bin/fits/fits.c b/bin/fits/fits.c
index e91adf9..97c8fd9 100644
--- a/bin/fits/fits.c
+++ b/bin/fits/fits.c
@@ -88,11 +88,11 @@ fits_print_extension_info(struct fitsparams *p)
 {
   uint16_t *ui16;
   fitsfile *fptr;
-  int hasblankname=0;
-  gal_data_t *cols=NULL, *tmp;
-  char **tstra, **estra, **sstra;
   size_t i, numext, *dsize, ndim;
+  int hascomments=0, hasblankname=0;
+  char **tstra, **estra, **sstra, **cmstra;
   int j, nc, numhdu, hdutype, status=0, type;
+  gal_data_t *tmp, *cols=NULL, *sizecol, *commentscol;
   char *msg, *tstr=NULL, sstr[1000], extname[FLEN_VALUE];
 
 
@@ -113,6 +113,8 @@ fits_print_extension_info(struct fitsparams *p)
   /* Allocate all the columns (in reverse order, since this is a simple
      linked list). */
   gal_list_data_add_alloc(&cols, NULL, GAL_TYPE_STRING, 1, &numext, NULL, 1,
+                          -1, 1, "HDU_COMMENT", "note", "Possible comment");
+  gal_list_data_add_alloc(&cols, NULL, GAL_TYPE_STRING, 1, &numext, NULL, 1,
                           -1, 1, "HDU_SIZE", "name", "Size of image or table "
                           "number of rows and columns.");
   gal_list_data_add_alloc(&cols, NULL, GAL_TYPE_STRING, 1, &numext, NULL, 1,
@@ -130,7 +132,12 @@ fits_print_extension_info(struct fitsparams *p)
   ui16  = cols->array;
   estra = cols->next->array;
   tstra = cols->next->next->array;
-  sstra = cols->next->next->next->array;
+
+  sizecol = cols->next->next->next;
+  sstra = sizecol->array;
+
+  commentscol = cols->next->next->next->next;
+  cmstra = commentscol->array;
 
   cols->next->disp_width=15;
   cols->next->next->disp_width=15;
@@ -182,6 +189,16 @@ fits_print_extension_info(struct fitsparams *p)
       status=0;
 
 
+      /* Check if its a healpix grid. */
+      if( gal_fits_hdu_is_healpix(fptr) )
+        {
+          hascomments=1;
+          gal_checkset_allocate_copy("HEALpix", cmstra+i);
+        }
+      else
+        gal_checkset_allocate_copy(GAL_BLANK_STRING, cmstra+i);
+
+
       /* Write the size into a string. 'sprintf' returns the number of
          written characters (excluding the '\0'). So for each dimension's
          size that is written, we add to 'nc' (the number of
@@ -229,6 +246,12 @@ fits_print_extension_info(struct fitsparams *p)
   /* Close the file. */
   fits_close_file(fptr, &status);
 
+  /* If there weren't any comments, don't print the comment column. */
+  if(hascomments==0)
+    {
+      sizecol->next=NULL;
+      gal_data_free(commentscol);
+    }
 
   /* Print the results. */
   if(!p->cp.quiet)
@@ -244,6 +267,9 @@ fits_print_extension_info(struct fitsparams *p)
       printf(" Column 3: Image data type or 'table' format (ASCII or "
              "binary).\n");
       printf(" Column 4: Size of data in HDU.\n");
+      if(hascomments)
+        printf(" Column 5: Comments about the HDU (e.g., if its HEALpix, or "
+               "etc).\n");
       printf("-----\n");
     }
   gal_table_write(cols, NULL, NULL, GAL_TABLE_FORMAT_TXT, NULL, NULL, 0);
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 8220749..93abd9d 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -8512,8 +8512,10 @@ $ astfits --write=MYKEY1,20.00,"An example keyword" 
--write=MYKEY2,fd
 @end example
 
 @cindex HDU
+@cindex HEALPix
 When no action is requested (and only a file name is given), Fits will print a 
list of information about the extension(s) in the file.
 This information includes the HDU number, HDU name (@code{EXTNAME} keyword), 
type of data (see @ref{Numeric data types}, and the number of data elements it 
contains (size along each dimension for images and table rows and columns).
+Optionally, a comment column is printed for special situations (like a 2D 
HEALPix grid that is usually stored as a 1D dataset/table).
 You can use this to get a general idea of the contents of the FITS file and 
what HDU to use for further processing, either with the Fits program or any 
other Gnuastro program.
 
 Here is one example of information about a FITS file with four extensions: the 
first extension has no data, it is a purely meta-data HDU (commonly used to 
keep meta-data about the whole file, or grouping of extensions, see @ref{Fits}).
@@ -23051,7 +23053,7 @@ Return the @code{DATASUM} of the given HDU in the given 
FITS file.
 For more on @code{DATASUM} in the FITS standard, see @ref{Keyword 
manipulation} (under the @code{checksum} component of @option{--write}).
 @end deftypefun
 
-@deftypefun {unsigned long} gal_fits_hdu_datasum_ptr (char @code{*fptr})
+@deftypefun {unsigned long} gal_fits_hdu_datasum_ptr (fitsfile @code{*fptr})
 @cindex @code{DATASUM}: FITS keyword
 Return the @code{DATASUM} of the already opened HDU in @code{fptr}.
 For more on @code{DATASUM} in the FITS standard, see @ref{Keyword 
manipulation} (under the @code{checksum} component of @option{--write}).
@@ -23062,6 +23064,12 @@ Return the format of the HDU as one of CFITSIO's 
recognized macros:
 @code{IMAGE_HDU}, @code{ASCII_TBL}, or @code{BINARY_TBL}.
 @end deftypefun
 
+@deftypefun int gal_fits_hdu_is_healpix (fitsfile @code{*fptr})
+@cindex HEALPix
+Return @code{1} if the dataset may be a HEALpix grid and @code{0} otherwise.
+Technically, it is considered to be a HEALPix if the HDU isn't an ASCII table, 
and has the @code{NSIDE}, @code{FIRSTPIX} and @code{LASTPIX}.
+@end deftypefun
+
 @deftypefun {fitsfile *} gal_fits_hdu_open (char @code{*filename}, char 
@code{*hdu}, int @code{iomode})
 Open the HDU/extension @code{hdu} from @file{filename} and return a pointer
 to CFITSIO's @code{fitsfile}. @code{iomode} determines how the FITS file
diff --git a/lib/fits.c b/lib/fits.c
index 023ecff..d894686 100644
--- a/lib/fits.c
+++ b/lib/fits.c
@@ -722,6 +722,30 @@ gal_fits_hdu_format(char *filename, char *hdu)
 
 
 
+/* Return 1 if the HDU appears to be a HEALpix grid. */
+int
+gal_fits_hdu_is_healpix(fitsfile *fptr)
+{
+  long value;
+  int hdutype, status=0;
+
+  /* An ASCII table can't be a healpix table. */
+  if (fits_get_hdu_type(fptr, &hdutype, &status) )
+    gal_fits_io_error(status, NULL);
+  if(hdutype==ASCII_TBL) return 0;
+
+  /* If all these keywords are present, then 'status' will be 0
+     afterwards. */
+  fits_read_key_lng(fptr, "NSIDE",    &value, 0x0, &status);
+  fits_read_key_lng(fptr, "FIRSTPIX", &value, 0x0, &status);
+  fits_read_key_lng(fptr, "LASTPIX",  &value, 0x0, &status);
+  return !status;
+}
+
+
+
+
+
 /* Open a given HDU and return the FITS pointer. 'iomode' determines how
    the FITS file will be opened: only to read or to read and write. You
    should use the macros given by the CFITSIO header:
@@ -837,8 +861,21 @@ gal_fits_hdu_open_format(char *filename, char *hdu, int 
img0_tab1)
   else
     {
       if(hdutype!=IMAGE_HDU)
-        error(EXIT_FAILURE, 0, "%s (hdu: %s): not an image",
-              filename, hdu);
+        {
+          /* Let the user know. */
+          if( gal_fits_hdu_is_healpix(fptr) )
+            error(EXIT_FAILURE, 0, "%s (hdu: %s): appears to be a HEALPix 
table "
+                  "(which is a 2D dataset on a spherical surface: stored as "
+                  "a 1D table). You can use the 'HPXcvt' command-line utility "
+                  "to convert it to a 2D image that can easily be used by "
+                  "other programs. 'HPXcvt' is built and installed as part of "
+                  "WCSLIB (which is a mandatory dependency of Gnuastro, so "
+                  "you should already have it), run 'man HPXcvt' for more",
+                  filename, hdu);
+          else
+            error(EXIT_FAILURE, 0, "%s (hdu: %s): not an image",
+                  filename, hdu);
+        }
     }
 
   /* Clean up and return. */
diff --git a/lib/gnuastro/fits.h b/lib/gnuastro/fits.h
index 9032db9..b8637b5 100644
--- a/lib/gnuastro/fits.h
+++ b/lib/gnuastro/fits.h
@@ -163,6 +163,9 @@ gal_fits_hdu_datasum_ptr(fitsfile *fptr);
 int
 gal_fits_hdu_format(char *filename, char *hdu);
 
+int
+gal_fits_hdu_is_healpix(fitsfile *fptr);
+
 fitsfile *
 gal_fits_hdu_open(char *filename, char *hdu, int iomode);