Smarter handling of map polarization properties

* Add TT, TQ, etc to MapPolType enumerator, and set the pol_type attribute for
  each map weights item.  Additionally, the TU and QU attributes of a
  G3SkyMapWeights object are multiplied by -1 if their polarization convention
  is changed from IAU to COSMO or vice versa (just like U maps).

* Where possible, determine parameters from the input map stub, similarly to the
  way these are handled in the MapBinner pipeline module.  Most noteably, the
  G3SkyMapWeights constructor now determines its polarization state based on the
  pol_conv property of the input reference map.  Similarly, the InjectMapStub
  pipeline module determines whether the output frame is weighted and/or
  polarized based on the properties of the input stub.

maps/include/maps/G3SkyMap.h

@@ -33,7 +33,18 @@ class G3SkyMap {
 		U = 2,
 		E = 3,
 		B = 4,
-		None = 7
+		None = 7,
+		TT = 8,
+		TQ = 9,
+		TU = 10,
+		QQ = 11,
+		QU = 12,
+		UU = 13,
+		TE = 14,
+		TB = 15,
+		EE = 16,
+		EB = 17,
+		BB = 18,
 	};
 
 	enum MapPolConv {
@@ -78,6 +89,8 @@ class G3SkyMap {
 	MapPolConv GetPolConv() const { return pol_conv_; };
 	void SetPolConv(MapPolConv pol_conv);
 
+	bool IsPolarized() const { return pol_conv_ != ConvNone; }
+
 	virtual bool IsCompatible(const G3SkyMap & other) const {
 		if (shape().size() != other.shape().size())
 			return false;
@@ -378,7 +391,7 @@ class G3SkyMapWeights : public G3FrameObject {
 	G3SkyMapWeights() {}
 
 	// Instantiate weight maps based on the metadata of a reference map
-	G3SkyMapWeights(G3SkyMapConstPtr ref_map, bool polarized = true);
+	G3SkyMapWeights(G3SkyMapConstPtr ref_map);
 
 	G3SkyMapWeights(const G3SkyMapWeights &r, bool copy_data = true);
 	G3SkyMapPtr TT, TQ, TU, QQ, QU, UU;

maps/python/fitsio.py

@@ -167,12 +167,11 @@ def load_skymap_fits(filename, hdu=None, keys=None, memmap=False, apply_units=Fa
             if map_type == 'flat' and hdr.get('ISWEIGHT', None):
                 # flat map weights
                 assert('T' in output)
-                if pol is None:
-                    pol = True if ('Q' in output and 'U' in output) else False
-                weight_map = output.setdefault(
-                    'W', G3SkyMapWeights(output['T'], polarized=pol)
-                )
-                fm = FlatSkyMap(data, **map_opts)
+                if 'W' not in output:
+                    output['W'] = G3SkyMapWeights()
+                weight_map = output['W']
+                pol_type = getattr(MapPolType, hdr['WTYPE'], None)
+                fm = FlatSkyMap(data, pol_type=pol_type, **map_opts)
                 fm.overflow = overflow
                 setattr(weight_map, hdr['WTYPE'], fm)
                 del data
@@ -242,24 +241,18 @@ def load_skymap_fits(filename, hdu=None, keys=None, memmap=False, apply_units=Fa
                         'TISWGT{:d}'.format(cidx + 1),
                         col in ['TT', 'TQ', 'TU', 'QQ', 'QU', 'UU'],
                     )
+                    pol_type = getattr(MapPolType, col, None)
+                    mdata = (pix, data, nside) if pix is not None else data
+
                     if weighted:
                         assert('T' in output)
-                        if pol is None:
-                            pol = True if ('Q' in output and 'U' in output) else False
-                        weight_map = output.setdefault(
-                            'W', G3SkyMapWeights(output['T'], polarized=pol)
-                        )
-
-                        mdata = (pix, data, nside) if pix is not None else data
-                        hm = HealpixSkyMap(mdata, **map_opts)
+                        if 'W' not in output:
+                            output['W'] = G3SkyMapWeights()
+                        weight_map = output['W']
+                        hm = HealpixSkyMap(mdata, pol_type=pol_type, **map_opts)
                         hm.overflow = overflow
-
                         setattr(weight_map, col, hm)
-
                     else:
-                        pol_type = getattr(MapPolType, col, None)
-                        mdata = (pix, data, nside) if pix is not None else data
-
                         hm = HealpixSkyMap(mdata, pol_type=pol_type, **map_opts)
                         output[col] = hm
 

maps/python/map_modules.py

@@ -157,21 +157,26 @@ def MakeMapsPolarized(frame, pol_conv=maps.MapPolConv.IAU):
 
     qmap = frame["T"].clone(False)
     qmap.pol_type = maps.MapPolType.Q
+    qmap.pol_conv = pol_conv
     frame["Q"] = qmap
     umap = frame["T"].clone(False)
     umap.pol_type = maps.MapPolType.U
     umap.pol_conv = pol_conv
     frame["U"] = umap
     mask = wgt.to_mask().to_map()
 
-    wgt_out = maps.G3SkyMapWeights(frame["T"], polarized=True)
+    wgt_out = maps.G3SkyMapWeights()
     wgt_out.TT = wgt
     wgt_out.TQ = wgt.clone(False)
     wgt_out.TU = wgt.clone(False)
     wgt_out.QQ = mask
     wgt_out.QU = wgt.clone(False)
     wgt_out.UU = mask.clone(True)
 
+    for k in wgt_out.keys():
+        wgt_out[k].pol_type = getattr(maps.MapPolType, k)
+        wgt_out[k].pol_conv = pol_conv
+
     frame["Wpol"] = wgt_out
 
     return frame
@@ -190,7 +195,7 @@ def MakeMapsUnpolarized(frame):
     del frame["Q"]
     del frame["U"]
 
-    wgt_out = maps.G3SkyMapWeights(frame["T"], polarized=False)
+    wgt_out = maps.G3SkyMapWeights()
     wgt_out.TT = wgt
 
     frame["Wunpol"] = wgt_out
@@ -354,30 +359,19 @@ class InjectMapStub(object):
     map_id : string
         Id to assign to the new map frame
     map_stub : G3SkyMap instance
-        Map stub from which to clone the Stokes maps and weights.
-    polarized : bool
-        If True, add Q and U maps to stub frame, and ensure that weights are
-        polarized.  Otherwise, only a T map is created.
-    weighted : bool
-        If True, add weights to the stub frame.
-    pol_conv : MapPolConv instance
-        Polarization convention to use.
+        Map stub from which to clone the Stokes maps and weights.  If the
+        `weighted` attribute of the stub is True, the output frame will include
+        weights.  If the `pol_conv` attribute of the stub is not None, the
+        output frame will include Q and U maps (and polarized weights).
     """
 
-    def __init__(
-        self,
-        map_id,
-        map_stub,
-        polarized=True,
-        weighted=True,
-        pol_conv=maps.MapPolConv.IAU,
-    ):
+    def __init__(self, map_id, map_stub):
         self.map_frame = core.G3Frame(core.G3FrameType.Map)
         self.map_frame["Id"] = map_id
 
         map_stub = map_stub.clone(False)
-        map_stub.weighted = weighted
-        map_stub.pol_conv = pol_conv
+        weighted = map_stub.weighted
+        polarized = map_stub.polarized
 
         T = map_stub.clone(False)
         T.pol_type = maps.MapPolType.T
@@ -390,7 +384,7 @@ def __init__(
             U.pol_type = maps.MapPolType.U
             self.map_frame["U"] = U
         if weighted:
-            W = maps.G3SkyMapWeights(map_stub, polarized)
+            W = maps.G3SkyMapWeights(map_stub)
             self.map_frame["Wpol" if polarized else "Wunpol"] = W
 
     def __call__(self, frame):
@@ -773,6 +767,9 @@ def __call__(self, frame):
                 "Coordinate rotation of polarized maps is not implemented"
             )
 
+        if "Q" in frame and not self.stub.polarized:
+            self.stub.pol_conv = frame["Q"].pol_conv
+
         for key in ["T", "Q", "U", "Wpol", "Wunpol", "H"]:
 
             if key not in frame:
@@ -785,7 +782,7 @@ def __call__(self, frame):
                 maps.reproj_map(m, mnew, rebin=self.rebin, interp=self.interp)
 
             elif key in ["Wpol", "Wunpol"]:
-                mnew = maps.G3SkyMapWeights(self.stub, key == "Wpol")
+                mnew = maps.G3SkyMapWeights(self.stub)
                 for wkey in mnew.keys():
                     maps.reproj_map(
                         m[wkey], mnew[wkey], rebin=self.rebin, interp=self.interp

maps/src/G3SkyMap.cxx

@@ -90,13 +90,28 @@ G3SkyMapWeights::serialize(A &ar, unsigned v)
 G3_SERIALIZABLE_CODE(G3SkyMap);
 G3_SERIALIZABLE_CODE(G3SkyMapWeights);
 
+
+static bool
+isumap(const G3SkyMap &m)
+{
+  if (m.pol_type == G3SkyMap::U)
+    return true;
+  if (m.pol_type == G3SkyMap::TU)
+    return true;
+  if (m.pol_type == G3SkyMap::QU)
+    return true;
+
+  return false;
+}
+
+
 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)
 {
-	if (pol_type == U && 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.");
 }
@@ -111,20 +126,20 @@ G3SkyMap::ArrayClone(boost::python::object val) const
 }
 
 
-G3SkyMapWeights::G3SkyMapWeights(G3SkyMapConstPtr ref, bool polarized) :
-    TT(ref->Clone(false)), TQ(polarized ? ref->Clone(false) : NULL),
-    TU(polarized ? ref->Clone(false) : NULL),
-    QQ(polarized ? ref->Clone(false) : NULL),
-    QU(polarized ? ref->Clone(false) : NULL),
-    UU(polarized ? ref->Clone(false) : NULL)
+G3SkyMapWeights::G3SkyMapWeights(G3SkyMapConstPtr ref) :
+    TT(ref->Clone(false)), TQ(ref->IsPolarized() ? ref->Clone(false) : NULL),
+    TU(ref->IsPolarized() ? ref->Clone(false) : NULL),
+    QQ(ref->IsPolarized() ? ref->Clone(false) : NULL),
+    QU(ref->IsPolarized() ? ref->Clone(false) : NULL),
+    UU(ref->IsPolarized() ? ref->Clone(false) : NULL)
 {
-	TT->pol_type = G3SkyMap::None;
-	if (polarized){
-		TQ->pol_type = G3SkyMap::None;
-		TU->pol_type = G3SkyMap::None;
-		QQ->pol_type = G3SkyMap::None;
-		QU->pol_type = G3SkyMap::None;
-		UU->pol_type = G3SkyMap::None;
+	TT->pol_type = G3SkyMap::TT;
+	if (ref->IsPolarized()){
+		TQ->pol_type = G3SkyMap::TQ;
+		TU->pol_type = G3SkyMap::TU;
+		QQ->pol_type = G3SkyMap::QQ;
+		QU->pol_type = G3SkyMap::QU;
+		UU->pol_type = G3SkyMap::UU;
 	}
 }
 
@@ -481,11 +496,11 @@ skymap_setitem(G3SkyMap &skymap, ssize_t i, double val)
 void
 G3SkyMap::SetPolConv(G3SkyMap::MapPolConv pol_conv)
 {
-	if (pol_type == G3SkyMap::U && pol_conv == G3SkyMap::ConvNone)
+	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 (pol_type != G3SkyMap::U ||
+	if (!isumap(*this) ||
 	    pol_conv == G3SkyMap::ConvNone ||
 	    pol_conv_ == G3SkyMap::ConvNone) {
 		pol_conv_ = pol_conv;
@@ -1426,6 +1441,17 @@ PYBINDINGS("maps") {
 	    .value("E", G3SkyMap::E)
 	    .value("B", G3SkyMap::B)
 	    .value("none", G3SkyMap::None) // "None" is reserved in python
+	    .value("TT", G3SkyMap::TT)
+	    .value("TQ", G3SkyMap::TQ)
+	    .value("TU", G3SkyMap::TU)
+	    .value("QQ", G3SkyMap::QQ)
+	    .value("QU", G3SkyMap::QU)
+	    .value("UU", G3SkyMap::UU)
+	    .value("TE", G3SkyMap::TE)
+	    .value("TB", G3SkyMap::TB)
+	    .value("EE", G3SkyMap::EE)
+	    .value("EB", G3SkyMap::EB)
+	    .value("BB", G3SkyMap::BB)
 	;
 	enum_none_converter::from_python<G3SkyMap::MapPolType>();
 
@@ -1453,6 +1479,8 @@ PYBINDINGS("maps") {
 	      "(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.")
+	    .add_property("polarized", &G3SkyMap::IsPolarized,
+	      "True if the pol_conv property is set to IAU or COSMO, False otherwise.")
 	    .def_readwrite("units", &G3SkyMap::units,
 	      "Unit class (core.G3TimestreamUnits) of the map (e.g. "
 	      "core.G3TimestreamUnits.Tcmb). Within each unit class, further "
@@ -1674,9 +1702,10 @@ PYBINDINGS("maps") {
 	boost::python::implicitly_convertible<G3SkyMapPtr, G3SkyMapConstPtr>();
 
 	EXPORT_FRAMEOBJECT(G3SkyMapWeights, init<>(),
-	    "Polarized (Mueller matrix) or unpolarized (scalar) map pixel weights.")
-	    .def(bp::init<G3SkyMapConstPtr, bool>(
-	      (bp::arg("skymap"), bp::arg("polarized") = true)))
+	    "Polarized (Mueller matrix) or unpolarized (scalar) map pixel weights."
+	    "Weights are polarized if the pol_conv attribute of the reference map is set.")
+	    .def(bp::init<G3SkyMapConstPtr>(
+	      (bp::arg("skymap"))))
 	    .def_readwrite("TT",&G3SkyMapWeights::TT, "Mueller matrix component map")
 	    .def_readwrite("TQ",&G3SkyMapWeights::TQ, "Mueller matrix component map")
 	    .def_readwrite("TU",&G3SkyMapWeights::TU, "Mueller matrix component map")

maps/src/MapBinner.cxx

@@ -152,8 +152,7 @@ MapBinner::MapBinner(std::string output_map_id, const G3SkyMap &stub_map,
 	T_ = stub_map.Clone(false);
 	T_->pol_type = G3SkyMap::T;
 	if (store_weight_map)
-		map_weights_ = G3SkyMapWeightsPtr(new G3SkyMapWeights(T_,
-		    T_->GetPolConv() != G3SkyMap::ConvNone));
+		map_weights_ = G3SkyMapWeightsPtr(new G3SkyMapWeights(T_));
 
 	if (T_->GetPolConv() != G3SkyMap::ConvNone) {
 		Q_ = stub_map.Clone(false);
@@ -207,8 +206,7 @@ MapBinner::Process(G3FramePtr frame, std::deque<G3FramePtr> &out)
 
 		if (map_weights_) {
 			out_frame->Put((Q_) ? "Wpol" : "Wunpol", map_weights_);
-			map_weights_ = G3SkyMapWeightsPtr(new G3SkyMapWeights(
-			    T_, T_->GetPolConv() != G3SkyMap::ConvNone));
+			map_weights_ = G3SkyMapWeightsPtr(new G3SkyMapWeights(T_));
 		}
 
 		out_frame->Put("StartTime", G3TimePtr(new G3Time(start_)));

maps/src/SingleDetectorBoresightBinner.cxx

@@ -109,6 +109,7 @@ SingleDetectorBoresightBinner::SingleDetectorBoresightBinner(
 {
 	template_ = stub_map.Clone(false);
 	template_->pol_type = G3SkyMap::T;
+	template_->SetPolConv(G3SkyMap::ConvNone);
 }
 
 void
@@ -172,7 +173,7 @@ SingleDetectorBoresightBinner::Process(G3FramePtr frame,
 			#endif
 		}
 		map_weights_ = G3SkyMapWeightsPtr(
-		    new G3SkyMapWeights(template_, false));
+		    new G3SkyMapWeights(template_));
 	} else {
 		if (template_->units != timestreams->GetUnits())
 			log_fatal("Timestreams have units that do not match "

maps/src/SingleDetectorMapBinner.cxx

@@ -57,6 +57,7 @@ SingleDetectorMapBinner::SingleDetectorMapBinner(
 {
 	template_ = stub_map.Clone(false);
 	template_->pol_type = G3SkyMap::T;
+	template_->SetPolConv(G3SkyMap::ConvNone);
 }
 
 void
@@ -120,7 +121,7 @@ SingleDetectorMapBinner::Process(G3FramePtr frame,
 
 			maps_[i.first].first = template_->Clone(false);
 			maps_[i.first].second = G3SkyMapWeightsPtr(
-			    new G3SkyMapWeights(template_, false));
+			    new G3SkyMapWeights(template_));
 #ifdef OPENMP_FOUND
 			// Create a list of detectors to satisfy OpenMP's need
 			// for scalar iteration.

maps/tests/map_modules.py

@@ -11,7 +11,7 @@
     pipe = core.G3Pipeline()
 
     pipe.Add(core.G3InfiniteSource, type=core.G3FrameType.Observation, n=1)
-    pipe.Add(maps.InjectMapStub, map_id="test_map", map_stub=m, polarized=False, weighted=True)
+    pipe.Add(maps.InjectMapStub, map_id="test_map", map_stub=m)
 
     def RandomMap(frame):
         if frame.type != core.G3FrameType.Map:

maps/tests/weights_operators.py

@@ -7,15 +7,16 @@
 
 for pol in [True, False]:
     # allocation
-    m = FlatSkyMap(500, 20, core.G3Units.arcmin, pol_conv=MapPolConv.IAU, proj=MapProjection.Proj5)
+    pol_conv = MapPolConv.IAU if pol else MapPolConv.none
+    m = FlatSkyMap(500, 20, core.G3Units.arcmin, pol_conv=pol_conv, proj=MapProjection.Proj5)
     mt = m.clone(False)
     mt.pol_type = MapPolType.T
     if pol:
         mq = m.clone(False)
         mq.pol_type = MapPolType.Q
         mu = m.clone(False)
         mu.pol_type = MapPolType.U
-    mw = G3SkyMapWeights(m, polarized=pol)
+    mw = G3SkyMapWeights(m)
     assert(mw.size == m.size)
     assert(mw.shape == m.shape)
     assert(mw.npix_allocated == 0)