diff --git a/maps/include/maps/G3SkyMap.h b/maps/include/maps/G3SkyMap.h index b9af0042..59673d99 100644 --- a/maps/include/maps/G3SkyMap.h +++ b/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; diff --git a/maps/python/fitsio.py b/maps/python/fitsio.py index 2b200573..6edee525 100644 --- a/maps/python/fitsio.py +++ b/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 diff --git a/maps/python/map_modules.py b/maps/python/map_modules.py index e91d40a1..38800665 100644 --- a/maps/python/map_modules.py +++ b/maps/python/map_modules.py @@ -157,6 +157,7 @@ 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 @@ -164,7 +165,7 @@ def MakeMapsPolarized(frame, pol_conv=maps.MapPolConv.IAU): 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) @@ -172,6 +173,10 @@ def MakeMapsPolarized(frame, pol_conv=maps.MapPolConv.IAU): 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 diff --git a/maps/src/G3SkyMap.cxx b/maps/src/G3SkyMap.cxx index a48b3620..9cb828be 100644 --- a/maps/src/G3SkyMap.cxx +++ b/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(); @@ -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(); EXPORT_FRAMEOBJECT(G3SkyMapWeights, init<>(), - "Polarized (Mueller matrix) or unpolarized (scalar) map pixel weights.") - .def(bp::init( - (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( + (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") diff --git a/maps/src/MapBinner.cxx b/maps/src/MapBinner.cxx index 9d5f31b6..386d951c 100644 --- a/maps/src/MapBinner.cxx +++ b/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 &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_))); diff --git a/maps/src/SingleDetectorBoresightBinner.cxx b/maps/src/SingleDetectorBoresightBinner.cxx index 2dcb590d..4f455d1a 100644 --- a/maps/src/SingleDetectorBoresightBinner.cxx +++ b/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 " diff --git a/maps/src/SingleDetectorMapBinner.cxx b/maps/src/SingleDetectorMapBinner.cxx index fbc7e848..2edfa881 100644 --- a/maps/src/SingleDetectorMapBinner.cxx +++ b/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. diff --git a/maps/tests/map_modules.py b/maps/tests/map_modules.py index de046437..827bd9ba 100644 --- a/maps/tests/map_modules.py +++ b/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: diff --git a/maps/tests/weights_operators.py b/maps/tests/weights_operators.py index 1886493a..d5268fb9 100755 --- a/maps/tests/weights_operators.py +++ b/maps/tests/weights_operators.py @@ -7,7 +7,8 @@ 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: @@ -15,7 +16,7 @@ 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)