utils/reproject.py: Add more logging and some minor changes

1 files changed, 22 insertions, 1 deletions

diff --git a/fg21sim/utils/reproject.py b/fg21sim/utils/reproject.py
index b129dc6..08dd1c9 100644
--- a/fg21sim/utils/reproject.py
+++ b/fg21sim/utils/reproject.py
@@ -3,6 +3,10 @@
 
 """
 FITS WCS (world coordinate system) reprojection utilities.
+
+zea2healpix:
+    Reproject the maps in ZEA (zenithal/azimuthal equal area) projection
+    to Galactic frame and organize in HEALPix format.
 """
 
 import logging
@@ -60,6 +64,15 @@ def _convert_wcs(lon_in, lat_in, frame_in, frame_out):
     else:
         coordframe_out, lon_out_unit, lat_out_unit = frame_out
     #
+    logger.info("Convert coordinates from {0} to {1}".format(coordframe_in,
+                                                             coordframe_out))
+    logger.info("Input coordinates units: "
+                "{0} (longitude), {1} (latitude)".format(lon_in_unit,
+                                                         lat_in_unit))
+    logger.info("Output coordinates units: "
+                "{0} (longitude), {1} (latitude)".format(lon_out_unit,
+                                                         lat_out_unit))
+    #
     data = UnitSphericalRepresentation(lon_in*lon_in_unit, lat_in*lat_in_unit)
     coords_in = coordframe_in.realize_frame(data)
     coords_out = coords_in.transform_to(coordframe_out)
@@ -117,29 +130,37 @@ def _image_to_healpix(image, wcs, nside, order=1, hemisphere=None):
     #
     npix = hp.nside2npix(nside)
     hpidx = np.arange(npix).astype(np.int)
+    logger.info("Output HEALPix: Nside={0}, Npixel={1}".format(nside, npix))
     # Calculate the longitude and latitude in frame of output HEALPix
+    logger.info("Calculate the longitudes and latitudes on the HEALPix grid")
     theta, phi = hp.pix2ang(nside, hpidx, nest=False)
     lon_hp = np.degrees(phi)
     lat_hp = 90.0 - np.degrees(theta)
     # Convert between the celestial coordinate systems
     coordsys_hp = Galactic()
+    logger.info("Output HEALPix uses frame: {0}".format(coordsys_hp))
     frame_hp = (coordsys_hp, au.deg, au.deg)
     lon_in, lat_in = _convert_wcs(lon_hp, lat_hp, frame_hp, wcs)
     # Filter the pixels on the specified hemisphere
     if hemisphere is None:
         mask = np.ones(lat_in.shape).astype(np.bool)
+        logger.info("NO hemisphere constraint specified")
     elif hemisphere.upper() in ["N", "NORTH", "NORTHERN"]:
         # northern hemisphere (include the equator)
         mask = lat_in >= 0.0
+        logger.info("Only process the NORTHERN hemisphere (include EQUATOR")
     else:
         # southern hemisphere
         mask = lat_in < 0.0
+        logger.info("Only process the SOUTHERN hemisphere (EXCLUDE equator")
     lon_in = lon_in[mask]
     lat_in = lat_in[mask]
     # Look up pixels in the input coordinate system
-    # NOTE: note the order of returns: (Y, X)
+    # XXX/NOTE: note the order of returns: (Y, X)
     yi, xi = wcs.wcs_world2pix(lon_in, lat_in, 0)
     # Interpolate to obtain the HEALPix data from the input image
+    logger.info("Calculate the HEALPix data by interpolating on input image")
+    logger.info("Interpolation order: {0}".format(order))
     data = map_coordinates(image, [xi, yi], order=order,
                            mode="constant", cval=np.nan)
     # Make the HEALPix array with above hemisphere mask considered