Merge branch 'polplots' into 'devel'

Polplots

See merge request parras/resolve!23

resolve/ubik_tools/plot_polarization.py

@@ -139,7 +139,7 @@ def _extent(sdom, offset=None):
         return [-xlim+ox, xlim+ox, -ylim+oy, ylim+oy]
 
 
-def polarization_quiver(ax, sky_field, nquivers=100, pfrac=0.001):
+def polarization_quiver(ax, sky_field, nquivers=100, pfrac=0.001, vmin=None, vmax=None):
     assert_sky_domain(sky_field.domain)
     pdom, tdom, fdom, sdom = sky_field.domain
     assert all((pol in pdom.labels) for pol in ["I", "Q", "U"])
@@ -147,7 +147,7 @@ def polarization_quiver(ax, sky_field, nquivers=100, pfrac=0.001):
     ang = polarization_angle(sky_field).val[0, 0]
     lin = linear_polarization(sky_field).val[0, 0]
 
-    im = ax.imshow(lin.T, cmap="inferno", norm=LogNorm(), origin="lower",
+    im = ax.imshow(lin.T, cmap="inferno", norm=LogNorm(vmin=vmin, vmax=vmax), origin="lower",
                    extent=_extent(sdom))
     scale = np.max(lin)*max(lin.shape) * 5
 
@@ -169,6 +169,24 @@ def polarization_quiver(ax, sky_field, nquivers=100, pfrac=0.001):
               scale=1.1/dst/skip, scale_units="xy",
               alpha=.6,
               )
+    plt.colorbar(im, ax=ax, orientation="horizontal", label="Linear polarized flux [Jy/sr]")
+
+
+def polarized_fraction(ax, sky_field, pfrac=0.001, vmax=None):
+    assert_sky_domain(sky_field.domain)
+    pdom, tdom, fdom, sdom = sky_field.domain
+    assert all((pol in pdom.labels) for pol in ["I", "Q", "U"])
+    assert tdom.size == fdom.size == 1
+    lin = linear_polarization(sky_field).val[0, 0]
+    total_int = sky_field.val[pdom.label2index("I"), 0, 0]
+
+    frac = lin/total_int
+    frac = np.ma.masked_where(lin < pfrac * np.max(lin), frac)
+
+    im = ax.imshow(frac.T, cmap="viridis_r", origin="lower", vmin=0, vmax=vmax,
+                   extent=_extent(sdom))
+
+    plt.colorbar(im, ax=ax, orientation="horizontal", label="Polarized fraction [1]")
 
 
 def polarization_angle(sky_field, faradaycorrection=0):

resolve/weighting_model.py

@@ -77,11 +77,10 @@ def log_weighting_model(cfg, obs, sky_domain):
                     tmpop.append(foo.ducktape(key).ducktape_left(key))
             linear_interpolation = reduce(add, tmpop)
             restructure = _CustomRestructure(linear_interpolation.target, oo.vis.domain)
-            ift.extra.check_linear_operator(restructure)
             tmpop = (restructure @ linear_interpolation @ log_weights).scale(-2)
             if oo.is_single_precision():
                 tmpop = DtypeConverter(tmpop.target, np.float64, np.float32) @ tmpop
-            tmpop = ift.Adder(oo.weight.log()) @ tmpop
+            tmpop = ift.Adder(_log_if_not_zero(oo.weight)) @ tmpop
             op.append(tmpop)
         return op, additional
     if cfg["model"] == "independent gamma":
@@ -126,3 +125,9 @@ class _CustomRestructure(ift.LinearOperator):
                 for ii in range(self.target.shape[2]):
                     res[_polfreq_key(pp, ii)] = x[pp, :, ii]
         return ift.makeField(self._tgt(mode), res)
+
+
+def _log_if_not_zero(fld):
+    res = fld.log().val_rw()
+    res[fld.val == 0.] = 0.
+    return ift.makeField(fld.domain, res)