Implement HpxNDMap.stack()
@@ -346,6 +346,36 @@ class HpxNDMap(HpxMap):
geom=geom, data=data, unit=self.unit, meta=self.meta
)
+ def stack(self, other, weights=None):
+ """Stack cutout into map.
+ Parameters
+ ----------
+ other : `HpxNDMap`
+ Other map to stack
+ weights : `HpxNDMap`
+ Array to be used as weights. The spatial geometry must be equivalent
+ to `other` and additional axes must be broadcastable.
+ """
+ if self.geom == other.geom:
+ idx = None
+ elif self.geom.is_aligned(other.geom):
+ if self.geom.is_allsky:
+ idx = other.geom._ipix
+ else:
+ idx = self.geom.to_image().global_to_local((other.geom._ipix,))[0]
+ else:
+ raise ValueError(
+ "Can only stack equivalent maps or cutout of the same map."
+ )
+
+ data = other.quantity.to_value(self.unit)
+
+ if weights is not None:
+ if not other.geom.to_image() == weights.geom.to_image():
+ raise ValueError("Incompatible spatial geoms between map and weights")
+ data = data * weights.data
+ self.data[..., idx] += data
+
def get_by_idx(self, idx):
# inherited docstring
idx = pix_tuple_to_idx(idx)
