New SetPolConv module for changing polarization convention

The G3SkyMap.pol_conv attribute is now a simple public attribute, and will not
result in changing the sign of the U map when updated.

maps/README.rst

@@ -27,7 +27,7 @@ The following attributes are common to all G3SkyMap subclasses:
   The Stokes polarization of the map object, which is an instance of the ``MapPolType`` enum, and can have the value ``T``, ``Q``, ``U`` or None.
 
 ``pol_conv``
-  The polarization convention used to encode the Q and U Stokes orientations relative to the coordinate axes.  This attribute is an instance of the ``MapPolConv`` enum, which can have the value ``IAU``, ``COSMO`` or None.  Both IAU and COSMO polarization conventions are supported in polarization-aware functions (e.g. ``FlattenPol``), but most default to using the IAU convention.  Warnings will be raised when a polarized map is used without a polarization convention set.  Changing the polarization convention between IAU and COSMO on a ``U`` map results in flipping the sign of all pixels in the map.
+  The polarization convention used to encode the Q and U Stokes orientations relative to the coordinate axes.  This attribute is an instance of the ``MapPolConv`` enum, which can have the value ``IAU``, ``COSMO`` or None.  Both IAU and COSMO polarization conventions are supported in polarization-aware functions (e.g. ``FlattenPol``), but most default to using the IAU convention.  Warnings will be raised when a polarized map is used without a polarization convention set.  Use the ``SetPolConv`` pipeline module to change the polarization convention between IAU and COSMO for an entire map frame.  This will result in flipping the sign of all pixels in the ``U`` map as well as the ``TU`` and ``QU`` weights.
 
 ``units``
   The units system in which the map is computed, stored as an instance of the ``G3TimestreamUnits`` enum, typically ``Tcmb``.

maps/include/maps/G3SkyMap.h

@@ -66,6 +66,7 @@ class G3SkyMap {
 	MapCoordReference coord_ref;
 	G3Timestream::TimestreamUnits units;
 	MapPolType pol_type;
+	MapPolConv pol_conv;
 	bool weighted;
 	double overflow;
 
@@ -86,10 +87,7 @@ class G3SkyMap {
 	virtual size_t NpixNonZero(void) const = 0;  // nonzero
 	virtual std::vector<size_t> shape(void) const = 0;  // map shape
 
-	MapPolConv GetPolConv() const { return pol_conv_; };
-	void SetPolConv(MapPolConv pol_conv);
-
-	bool IsPolarized() const { return pol_conv_ != ConvNone; }
+	bool IsPolarized() const { return pol_conv != ConvNone; }
 
 	virtual bool IsCompatible(const G3SkyMap & other) const {
 		if (shape().size() != other.shape().size())
@@ -218,8 +216,6 @@ class G3SkyMap {
 	    bool zero_infs = false) const;
 
 protected:
-	MapPolConv pol_conv_;
-
 	virtual void InitFromV1Data(std::vector<size_t>,
 	    const std::vector<double> &) {
 		throw std::runtime_error("Initializing from V1 not implemented");

maps/python/map_modules.py

@@ -5,6 +5,7 @@
     "CompactMaps",
     "RemoveWeights",
     "ApplyWeights",
+    "SetPolConv",
     "FlattenPol",
     "MakeMapsPolarized",
     "MakeMapsUnpolarized",
@@ -101,6 +102,53 @@ def ApplyWeights(frame):
     return frame
 
 
+@core.indexmod
+def SetPolConv(frame, pol_conv=maps.MapPolConv.IAU):
+    """
+    Set or change the polarization convention of the input polarized map frame.
+    If switching between IAU and COSMO conventions, flip the sign of the U map
+    and the TU and QU weights.  Otherwise, just set the polarization convention
+    for all maps and weights.
+    """
+    if not ("Q" in frame and "U" in frame):
+        # only polarized frames
+        return frame
+
+    if pol_conv == maps.MapPolConv.none or pol_conv is None:
+        raise ValueError("Polarized maps must have pol_conv set to IAU or COSMO")
+
+    tmap = frame.pop("T")
+    tmap.pol_conv = pol_conv
+
+    qmap = frame.pop("Q")
+    qmap.pol_conv = pol_conv
+
+    umap = frame.pop("U")
+    flip = umap.polarized and umap.pol_conv != pol_conv
+    umap.pol_conv = pol_conv
+
+    wmap = None
+    if "Wpol" in frame:
+        wmap = frame.pop("Wpol")
+        for k in wmap.keys():
+            wmap[k].pol_conv = pol_conv
+
+    # flip sign if switching conventions
+    if flip:
+        umap *= -1.0
+        if wmap is not None:
+            wmap.TU *= -1.0
+            wmap.QU *= -1.0
+
+    frame["T"] = tmap
+    frame["Q"] = qmap
+    frame["U"] = umap
+    if wmap is not None:
+        frame["Wpol"] = wmap
+
+    return frame
+
+
 @core.indexmod
 def FlattenPol(frame, invert=False):
     """
@@ -167,6 +215,8 @@ def MakeMapsPolarized(frame, pol_conv=maps.MapPolConv.IAU):
     umap.pol_conv = pol_conv
     frame["U"] = umap
     mask = wgt.to_mask().to_map()
+    mask.units = qmap.units
+    mask.weighted = qmap.weighted
 
     wgt_out = maps.G3SkyMapWeights()
     wgt_out.TT = wgt
@@ -248,6 +298,7 @@ def ValidateMaps(frame, ignore_missing_weights=False):
                 "Map frame %s: Map %s not compatible with T map" % (map_id, k),
                 unit="ValidateMaps",
             )
+
         if k in "TQU":
             if k == "U":
                 if isinstance(frame[k], maps.FlatSkyMap) and (
@@ -257,9 +308,10 @@ def ValidateMaps(frame, ignore_missing_weights=False):
                         "Map frame %s: Q and U maps have different flat_pol" % map_id,
                         unit="ValidateMaps",
                     )
-                if frame[k].pol_conv is maps.MapPolConv.none:
+            if k in "QU":
+                if not frame[k].polarized:
                     core.log_warn(
-                        "Map frame %s: U map polarization convention not set" % map_id,
+                        "Map frame %s: %s map polarization convention not set" % (map_id, k),
                         unit="ValidateMaps",
                     )
             if frame[k].weighted and not ignore_missing_weights:
@@ -280,28 +332,33 @@ def ValidateMaps(frame, ignore_missing_weights=False):
                             "Map frame %s: Missing polarized weights" % map_id,
                             unit="ValidateMaps",
                         )
-        else:
-            if frame[k].polarized and ("Q" not in frame or "U" not in frame):
-                core.log_fatal(
-                    "Map frame %s: Found unpolarized maps with polarized weights"
-                    % map_id,
-                    unit="ValidateMaps",
-                )
-            elif not frame[k].polarized and ("Q" in frame or "U" in frame):
+
+        elif k == "H":
+            continue
+
+        elif frame[k].polarized:
+            if "Q" not in frame or "U" not in frame:
                 core.log_fatal(
-                    "Map frame %s: Found polarized maps with unpolarized weights"
-                    % map_id,
+                    "Map frame %s: Found unpolarized maps with polarized weights" % map_id,
                     unit="ValidateMaps",
                 )
-            if (
-                frame[k].polarized
-                and isinstance(frame[k].QQ, maps.FlatSkyMap)
-                and frame[k].flat_pol != frame["Q"].flat_pol
-            ):
+            if frame[k].pol_conv != frame["U"].pol_conv:
                 core.log_fatal(
-                    "Map frame %s: %s and U maps have different flat_pol" % (map_id, k),
+                    "Map frame %s: %s and U maps have different pol_conv" % (map_id, k),
                     unit="ValidateMaps",
                 )
+            if isinstance(frame[k].QQ, maps.FlatSkyMap):
+                if frame[k].flat_pol != frame["Q"].flat_pol:
+                    core.log_fatal(
+                        "Map frame %s: %s and U maps have different flat_pol" % (map_id, k),
+                        unit="ValidateMaps",
+                    )
+
+        elif "Q" in frame or "U" in frame:
+            core.log_fatal(
+                "Map frame %s: Found polarized maps with unpolarized weights" % map_id,
+                unit="ValidateMaps",
+            )
 
 
 @core.indexmod

maps/python/skymapaddons.py

@@ -177,6 +177,8 @@ def skymapweights_reshape(self, width, height, fill=0):
 def skymapweights_getattr(self, x):
     if x == "flat_pol" and isinstance(self.TQ, FlatSkyMap):
         return getattr(self.TQ, x)
+    if x == "pol_conv" and isinstance(self.TU, FlatSkyMap):
+        return getattr(self.TU, x)
     if hasattr(self.TT, x):
         return getattr(self.TT, x)
     raise AttributeError("'{}' object has no attribute '{}'".format(type(self), x))

maps/src/FlatSkyMap.cxx

@@ -349,7 +349,7 @@ FlatSkyMap::Clone(bool copy_data) const
 		return boost::make_shared<FlatSkyMap>(*this);
 	else
 		return boost::make_shared<FlatSkyMap>(proj_info,
-		    coord_ref, weighted, units, pol_type, flat_pol_, pol_conv_);
+		    coord_ref, weighted, units, pol_type, flat_pol_, pol_conv);
 }
 
 double
@@ -567,7 +567,7 @@ FlatSkyMap::Description() const
 	default:
 		os << "unknown";
 	}
-	switch (pol_conv_) {
+	switch (pol_conv) {
 	case IAU:
 		os << " IAU";
 		break;
@@ -793,7 +793,7 @@ G3SkyMapPtr FlatSkyMap::Rebin(size_t scale, bool norm) const
 
 	FlatSkyProjection p(proj_info.Rebin(scale));
 	FlatSkyMapPtr out(new FlatSkyMap(p, coord_ref, weighted, units, pol_type,
-	    flat_pol_, pol_conv_));
+	    flat_pol_, pol_conv));
 
 	if (dense_)
 		out->ConvertToDense();
@@ -821,7 +821,7 @@ G3SkyMapPtr FlatSkyMap::ExtractPatch(size_t x0, size_t y0, size_t width, size_t
 
 	FlatSkyProjection p(proj_info.OverlayPatch(x0, y0, width, height));
 	FlatSkyMapPtr out(new FlatSkyMap(p, coord_ref, weighted, units, pol_type,
-	    flat_pol_, pol_conv_));
+	    flat_pol_, pol_conv));
 
 	if (fill != 0 && (width > xpix_ || height > ypix_))
 		(*out) += fill;

maps/src/G3SkyMap.cxx

@@ -48,9 +48,9 @@ G3SkyMap::serialize(A &ar, unsigned v)
 	ar & cereal::make_nvp("weighted", weighted);
 
 	if (v > 2) {
-		ar & cereal::make_nvp("pol_conv", pol_conv_);
+		ar & cereal::make_nvp("pol_conv", pol_conv);
 	} else {
-		pol_conv_ = ConvNone;
+		pol_conv = ConvNone;
 	}
 
 }
@@ -109,9 +109,9 @@ G3SkyMap::G3SkyMap(MapCoordReference coords, bool weighted,
     G3Timestream::TimestreamUnits units, MapPolType pol_type,
     MapPolConv pol_conv) :
     coord_ref(coords), units(units), pol_type(pol_type),
-    weighted(weighted), overflow(0), pol_conv_(pol_conv)
+    pol_conv(pol_conv), weighted(weighted), overflow(0)
 {
-	if (isumap(*this) && pol_conv_ == ConvNone)
+	if (isumap(*this) && pol_conv == ConvNone)
 		log_warn("Map object has pol_type U and unknown pol_conv. "
 			 "Set the pol_conv attribute to IAU or COSMO.");
 }
@@ -493,30 +493,6 @@ skymap_setitem(G3SkyMap &skymap, ssize_t i, double val)
 	skymap[i] = val;
 }
 
-void
-G3SkyMap::SetPolConv(G3SkyMap::MapPolConv pol_conv)
-{
-	if (isumap(*this) && pol_conv == G3SkyMap::ConvNone)
-		log_warn("Map object has pol_type U and unknown pol_conv. "
-			 "Set the pol_conv attribute to IAU or COSMO.");
-
-	if (!isumap(*this) ||
-	    pol_conv == G3SkyMap::ConvNone ||
-	    pol_conv_ == G3SkyMap::ConvNone) {
-		pol_conv_ = pol_conv;
-		return;
-	}
-
-	if (pol_conv != pol_conv_) {
-		log_warn("Sign of U map flipped in changing pol_conv from %s to %s",
-		    pol_conv == G3SkyMap::IAU ? "COSMO" : "IAU",
-		    pol_conv == G3SkyMap::IAU ? "IAU" : "COSMO");
-		(*this) *= -1;
-        }
-
-	pol_conv_ = pol_conv;
-}
-
 static bp::tuple
 skymap_shape(const G3SkyMap &skymap)
 {
@@ -1474,11 +1450,9 @@ PYBINDINGS("maps") {
 	    .def_readwrite("pol_type", &G3SkyMap::pol_type,
 	      "Polarization type (maps.MapPolType) of the map "
 	      "(e.g. maps.MapPolType.Q).")
-	    .add_property("pol_conv", &G3SkyMap::GetPolConv, &G3SkyMap::SetPolConv,
+	    .def_readwrite("pol_conv", &G3SkyMap::pol_conv,
 	      "Polarization convention (maps.MapPolConv) of the map "
-	      "(e.g. maps.MapPolConv.IAU or maps.MapPolConv.COSMO). "
-	      "Switching between IAU and COSMO conventions for a U map "
-	      "multiplies the U map by -1.")
+	      "(e.g. maps.MapPolConv.IAU or maps.MapPolConv.COSMO).")
 	    .add_property("polarized", &G3SkyMap::IsPolarized,
 	      "True if the pol_conv property is set to IAU or COSMO, False otherwise.")
 	    .def_readwrite("units", &G3SkyMap::units,

maps/src/G3SkyMapMask.cxx

@@ -15,7 +15,7 @@ G3SkyMapMask::G3SkyMapMask(const G3SkyMap &parent, bool use_data,
 	G3SkyMapPtr tmp = parent.Clone(false);
 	tmp->units = G3Timestream::None;
 	tmp->pol_type = G3SkyMap::None;
-	tmp->SetPolConv(G3SkyMap::ConvNone);
+	tmp->pol_conv = G3SkyMap::ConvNone;
 	tmp->weighted = false;
 	parent_ = tmp;
 	data_ = std::vector<bool>(parent.size());
@@ -31,7 +31,7 @@ G3SkyMapMask::G3SkyMapMask(const G3SkyMap &parent, boost::python::object v,
 	G3SkyMapPtr tmp = parent.Clone(false);
 	tmp->units = G3Timestream::None;
 	tmp->pol_type = G3SkyMap::None;
-	tmp->SetPolConv(G3SkyMap::ConvNone);
+	tmp->pol_conv = G3SkyMap::ConvNone;
 	tmp->weighted = false;
 	parent_ = tmp;
 	data_ = std::vector<bool>(parent.size());

maps/src/HealpixSkyMap.cxx

@@ -576,7 +576,7 @@ HealpixSkyMap::Clone(bool copy_data) const
 		return boost::make_shared<HealpixSkyMap>(*this);
 	else
 		return boost::make_shared<HealpixSkyMap>(nside(), weighted,
-		    nested(), coord_ref, units, pol_type, info_.shifted(), pol_conv_);
+		    nested(), coord_ref, units, pol_type, info_.shifted(), pol_conv);
 }
 
 double
@@ -869,7 +869,7 @@ HealpixSkyMap::Description() const
 	default:
 		os << "unknown";
 	}
-	switch (pol_conv_) {
+	switch (pol_conv) {
 	case IAU:
 		os << " IAU";
 		break;
@@ -1059,7 +1059,7 @@ G3SkyMapPtr HealpixSkyMap::Rebin(size_t scale, bool norm) const
 		return Clone(true);
 
 	HealpixSkyMapPtr out(new HealpixSkyMap(nside()/scale, weighted,
-	    nested(), coord_ref, units, pol_type, info_.shifted(), pol_conv_));
+	    nested(), coord_ref, units, pol_type, info_.shifted(), pol_conv));
 
 	if (dense_)
 		out->ConvertToDense();

maps/src/HitsBinner.cxx

@@ -117,7 +117,7 @@ HitsBinner::HitsBinner(std::string output_map_id, const G3SkyMap &stub_map,
 {
 	H_ = stub_map.Clone(false);
 	H_->pol_type = G3SkyMap::None;
-	H_->SetPolConv(G3SkyMap::ConvNone);
+	H_->pol_conv = G3SkyMap::ConvNone;
 	H_->units = G3Timestream::None;
 	H_->weighted = false;
 

maps/src/MapBinner.cxx

@@ -154,7 +154,7 @@ MapBinner::MapBinner(std::string output_map_id, const G3SkyMap &stub_map,
 	if (store_weight_map)
 		map_weights_ = G3SkyMapWeightsPtr(new G3SkyMapWeights(T_));
 
-	if (T_->GetPolConv() != G3SkyMap::ConvNone) {
+	if (T_->IsPolarized()) {
 		Q_ = stub_map.Clone(false);
 		Q_->pol_type = G3SkyMap::Q;
 		U_ = stub_map.Clone(false);

maps/src/MapMockObserver.cxx

@@ -112,9 +112,9 @@ MapMockObserver::MapMockObserver(std::string pointing, std::string timestreams,
 		log_fatal("If simulating polarized maps, pass both Q and U.");
 
 	if (U_) {
-		if (U_->GetPolConv() == G3SkyMap::ConvNone)
+		if (!(U_->IsPolarized()))
 			log_fatal("Missing pol_conv");
-		pol_sign_ = (U_->GetPolConv() == G3SkyMap::COSMO) ? -1 : 1;
+		pol_sign_ = (U_->pol_conv == G3SkyMap::COSMO) ? -1 : 1;
 	}
 }
 

maps/src/SingleDetectorBoresightBinner.cxx

@@ -109,7 +109,7 @@ SingleDetectorBoresightBinner::SingleDetectorBoresightBinner(
 {
 	template_ = stub_map.Clone(false);
 	template_->pol_type = G3SkyMap::T;
-	template_->SetPolConv(G3SkyMap::ConvNone);
+	template_->pol_conv = G3SkyMap::ConvNone;
 }
 
 void

maps/src/SingleDetectorMapBinner.cxx

@@ -57,7 +57,7 @@ SingleDetectorMapBinner::SingleDetectorMapBinner(
 {
 	template_ = stub_map.Clone(false);
 	template_->pol_type = G3SkyMap::T;
-	template_->SetPolConv(G3SkyMap::ConvNone);
+	template_->pol_conv = G3SkyMap::ConvNone;
 }
 
 void