From 1e6b0669374276a58cf1b081fa6efbdffd8266c2 Mon Sep 17 00:00:00 2001
From: Aaron LI <aaronly.me@outlook.com>
Date: Mon, 10 Oct 2016 17:12:48 +0800
Subject: utils: zea2healpix() also return the mask array

The mask array has values 0, 1, and 2, which indicate the different
statuses of the reprojected HEALPix pixels.

Also update the executable script to use the new return results.
---
 bin/zea2healpix            | 10 +++++-----
 fg21sim/utils/reproject.py | 18 +++++++++++++-----
 2 files changed, 18 insertions(+), 10 deletions(-)

diff --git a/bin/zea2healpix b/bin/zea2healpix
index c8879ab..1eb8730 100755
--- a/bin/zea2healpix
+++ b/bin/zea2healpix
@@ -74,11 +74,11 @@ def main():
         'See also {0}'.format(fg21sim.__url__)
     ]
 
-    hp_data, hp_header = zea2healpix(args.infile1, args.infile2,
-                                     nside=args.nside,
-                                     order=args.interp_order,
-                                     append_history=history,
-                                     append_comment=comments)
+    hp_data, hp_header, __ = zea2healpix(args.infile1, args.infile2,
+                                         nside=args.nside,
+                                         order=args.interp_order,
+                                         append_history=history,
+                                         append_comment=comments)
     write_fits_healpix(args.outfile, hpmap=hp_data, header=hp_header,
                        clobber=args.clobber)
     logger.info("HEALPix data write to FITS file: %s" % args.outfile)
diff --git a/fg21sim/utils/reproject.py b/fg21sim/utils/reproject.py
index 08dd1c9..f44c265 100644
--- a/fg21sim/utils/reproject.py
+++ b/fg21sim/utils/reproject.py
@@ -196,6 +196,14 @@ def zea2healpix(img1, img2, nside, order=1, inpaint=False,
         Reprojected HEALPix data
     hp_header : `~astropy.io.fits.Header`
         FITS header for the reprojected HEALPix data
+    hp_mask : 1D `~numpy.ndarray`
+        Array of same shape as the above `hp_data` indicating the status of
+        each pixel of the output array.
+        Values of "0" indicate the missing pixels (i.e., there is no
+        transformation to the input images); values of "1" indicate the output
+        pixel maps to one and only one of the input images; values of "2"
+        indicate the duplicate/overlapping pixels that map to both of the two
+        input images.
 
     NOTE
     ----
@@ -230,18 +238,18 @@ def zea2healpix(img1, img2, nside, order=1, inpaint=False,
     hp_data2 = _image_to_healpix(zea_img2, zea_wcs2, nside=nside, order=order,
                                  hemisphere=zea_hemisphere2.upper())
     # Merge the two HEALPix data
-    hp_footprint1 = (~np.isnan(hp_data1)).astype(np.int)
-    hp_footprint2 = (~np.isnan(hp_data2)).astype(np.int)
+    hp_mask = ((~np.isnan(hp_data1)).astype(np.int) +
+               (~np.isnan(hp_data2)).astype(np.int))
     hp_data1[np.isnan(hp_data1)] = 0.0
     hp_data2[np.isnan(hp_data2)] = 0.0
     hp_data = hp_data1 + hp_data2
     logger.info("Done reprojection and merge two hemispheres")
     # Duplicated pixels and missing pixels
-    pix_dup = (hp_footprint1 + hp_footprint2) == 2
+    pix_dup = (hp_mask == 2)
     if pix_dup.sum() > 0:
         logger.warning("Reprojected HEALPix data has %d duplicated pixel(s)" %
                        pix_dup.sum())
-    pix_missing = (hp_footprint1 + hp_footprint2) == 0
+    pix_missing = (hp_mask == 0)
     if pix_missing.sum() > 0:
         logger.warning("Reprojected HEALPix data has %d missing pixel(s)" %
                        pix_missing.sum())
@@ -258,4 +266,4 @@ def zea2healpix(img1, img2, nside, order=1, inpaint=False,
                                      append_history=append_history,
                                      append_comment=append_comment)
     logger.info("Made HEALPix FITS header")
-    return (hp_data, hp_header)
+    return (hp_data, hp_header, hp_mask)
-- 
cgit v1.2.2