[gnuastro-commits] master 19977d4: Fits: pixel scales also shown in arcsec/pixel

[Mohammad Akhlaghi]

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

    Fits: pixel scales also shown in arcsec/pixel
    
    Until now, the recently added '--pixelscale' option would only print the
    pixel scale in units of 'deg/pixel', which is the raw units of the CDELT
    keyword in the FITS WCS standrad.
    
    However, when the WCS axis is in units of degrees, it is most common to
    express the pixel scale in units of 'arcsec/pixel'. So besides printing the
    'deg/pixel', when the units of the axis is in units of degrees, the pixel
    scale is also printed in units of 'arcsec/pixel'.
    
    Also, with this commit, the pixel area (on an image, or on each slice of a
    3D cube) is printed and when we have a 3D dataset, the pixel volume is also
    printed. Similar to above, when the units are in 'deg', values are also
    printed in 'arcsec' and when we have meters ('m'), the value is also
    printed in Angestroms.
---
 NEWS              |  5 ++++
 bin/fits/fits.c   | 71 ++++++++++++++++++++++++++++++++++++++++++++++++++-----
 doc/gnuastro.texi | 12 ++++++++--
 3 files changed, 80 insertions(+), 8 deletions(-)

diff --git a/NEWS b/NEWS
index 50a159d..29c4912 100644
--- a/NEWS
+++ b/NEWS
@@ -99,6 +99,11 @@ See the end of the file for license conditions.
      added "Memory management" section of the book for a complete
      explanation of Gnuastro's new memory management strategy.
 
+  Fits:
+   - 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.
+
   Table:
    - Column arithmetic operators 'degree-to-ra' and 'degree-to-dec' will
      return the sexagesimal format of '_h_m_s' and '_d_m_s'
diff --git a/bin/fits/fits.c b/bin/fits/fits.c
index 90a9348..2a6e386 100644
--- a/bin/fits/fits.c
+++ b/bin/fits/fits.c
@@ -308,8 +308,8 @@ fits_pixelscale(struct fitsparams *p)
 {
   int nwcs=0;
   size_t i, ndim=0;
-  double *pixelscale;
   struct wcsprm *wcs;
+  double multip, *pixelscale;
 
   /* Read the desired WCS. */
   wcs=gal_wcs_read(p->filename, p->cp.hdu, 0, 0, &nwcs);
@@ -334,15 +334,74 @@ fits_pixelscale(struct fitsparams *p)
              "with '--quiet' or '-q')\n");
       printf("  Input: %s (hdu %s) has %zu dimensions.\n", p->filename,
              p->cp.hdu, ndim);
-      printf("  Pixel scale in each dimension:\n");
+      printf("  Pixel scale in each FITS dimension:\n");
       for(i=0;i<ndim;++i)
-        printf("    %zu: %g (%s/pixel)\n", i+1, pixelscale[i], wcs->cunit[i]);
+        {
+          if( !strcmp(wcs->cunit[i], "deg") )
+            printf("    %zu: %g (%s/pixel) = %g (arcsec/pixel)\n", i+1,
+                   pixelscale[i], wcs->cunit[i], pixelscale[i]*3600);
+          else
+            printf("    %zu: %g (%s/slice)\n", i+1,
+                   pixelscale[i], wcs->cunit[i]);
+        }
+
+      /* Pixel area/volume. */
+      if(ndim>=2)
+        {
+          /* Multiply the values in each dimension. */
+          multip=pixelscale[0]*pixelscale[1];
+
+          /* Pixel scale (applicable to 2 or 3 dimensions). */
+          printf("  Pixel area%s:\n", ndim==2?"":" (on each 2D slice) ");
+          if( strcmp(wcs->cunit[0], wcs->cunit[1]) )
+            printf("    %g (%s*%s)\n", multip, wcs->cunit[0],
+                   wcs->cunit[1]);
+          else
+            if( strcmp(wcs->cunit[0], "deg") )
+              printf("    %g (%s^2)\n", multip, wcs->cunit[0]);
+            else
+              printf("    %g (deg^2) = %g (arcsec^2)\n",
+                     multip, multip*3600*3600 );
+
+          /* For a 3 dimensional dataset, we need need extra info. */
+          if(ndim>=3)
+            {
+              multip*=pixelscale[2];
+              printf("  Voxel volume:\n");
+              if( strcmp(wcs->cunit[0], wcs->cunit[1]) )
+                printf("    %g (%s*%s*%s)\n", multip, wcs->cunit[0],
+                       wcs->cunit[1], wcs->cunit[2]);
+              else
+                if( strcmp(wcs->cunit[0], "deg") )
+                  printf("    %g (%s^2*%s)\n", multip, wcs->cunit[0],
+                         wcs->cunit[2]);
+                else
+                  {
+                    if( strcmp(wcs->cunit[2], "m") )
+                      printf("    %g (deg^2*%s) = %g (arcsec^2*%s)\n",
+                             multip, wcs->cunit[2], multip*3600*3600,
+                             wcs->cunit[2]);
+                    else
+                      printf("    %g (deg^2*m) = %g (arcsec^2*m) = "
+                             "%g (arcsec^2*A)\n", multip, multip*3600*3600,
+                             multip*3600*3600*1e10);
+                  }
+            }
+        }
     }
   else
     {
-      for(i=0;i<ndim-1;++i)
-        printf("%g ", pixelscale[i]);
-      printf("%g\n", pixelscale[ndim-1]);
+      multip=1;
+      for(i=0;i<ndim;++i)
+        {
+          multip*=pixelscale[i];
+          printf("%g ", pixelscale[i]);
+        }
+      switch(ndim)
+        {
+        case 2: printf("%g\n", multip); break;
+        case 3: printf("%g %g\n", pixelscale[0]*pixelscale[1], multip); break;
+        }
     }
 
   /* Clean up. */
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 67d5222..b340996 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -8456,9 +8456,17 @@ For a complete list of basic meta-data on the extensions 
in a FITS file, don't u
 For more, see @ref{Invoking astfits}.
 
 @item --pixelscale
-Print the HDU's pixel-scale (change in world coordinate for one pixel along 
each dimension).
-Without the @option{--quiet} option, the output of @option{--pixelscale} is 
more human-friendly by printing the file/HDU name, number of dimensions, and 
the units of each number along with the actual pixel scales.
+Print the HDU's pixel-scale (change in world coordinate for one pixel along 
each dimension) and pixel area or voxel volume.
+Without the @option{--quiet} option, the output of @option{--pixelscale} has 
multiple lines and explanations, thus being more human-friendly.
+It prints the file/HDU name, number of dimensions, and the units along with 
the actual pixel scales.
+Also, when any of the units are in degrees, the pixel scales and area/volume 
are also printed in units of arcsecs.
+For 3D datasets, the pixel area (on each 2D slice of the 3D cube) is printed 
as well as the voxel volume.
+
 However, in scripts (that are to be run automatically), this human-friendly 
format is annoying, so when called with the @option{--quiet} option, only the 
pixel-scale value(s) along each dimension is(are) printed in one line.
+These numbers are followed by the pixel area (in the raw WCS units).
+For 3D datasets, this will be area on each 2D slice.
+Finally, for 3D datasets, a final number (the voxel volume) is printed.
+As a summary, in @option{--quiet} mode, for 2D datasets three numbers are 
printed and for 3D datasets, 5 numbers are printed.
 
 @item --skycoverage
 @cindex Image's sky coverage