Make software compatible with Visibilities class (#67)

* Fixed bug in the default creation of VisMeta * Made the software compatible with new Visibilities class --------- Co-authored-by: Shane Maloney <[email protected]>
TCDSolar · Jun 1, 2024 · d195e70 · d195e70
1 parent dcd2a3e
commit d195e70
Show file tree

Hide file tree

Showing 13 changed files with 225 additions and 210 deletions.
diff --git a/changelog/67.feature.rst b/changelog/67.feature.rst
@@ -0,0 +1 @@
+Make the software compatible with :class:`~xrayvision.visibility.Visibilities`.
diff --git a/examples/rhessi.py b/examples/rhessi.py
@@ -13,8 +13,8 @@
 
 from xrayvision.clean import vis_clean
 from xrayvision.imaging import vis_psf_map, vis_to_map
-from xrayvision.mem import mem
-from xrayvision.visibility import Visibility
+from xrayvision.mem import mem, resistant_mean
+from xrayvision.visibility import Visibilities, VisMeta
 
 ###############################################################################
 # We will use `astropy.io.fits` to download and open the RHESSI visibility fits
@@ -50,15 +50,17 @@
 ###############################################################################
 # Now we can create the visibility object from the filtered visibilities.
 
+meta = VisMeta({"vis_labels": vis_data["isc"]})
+
 vunit = apu.Unit("photon/(cm**2 s)")
-vis = Visibility(
-    vis=vis_data["obsvis"] * vunit,
+vis = Visibilities(
+    visibilities=vis_data["obsvis"] * vunit,
     u=vis_data["u"] / apu.arcsec,
     v=vis_data["v"] / apu.arcsec,
-    offset=vis_data["xyoffset"][0] * apu.arcsec,
+    phase_center=vis_data["xyoffset"][0] * apu.arcsec,
+    meta=meta,
+    amplitude_uncertainty=vis_data["sigamp"] * vunit,
 )
-setattr(vis, "amplitude_error", vis_data["sigamp"] * vunit)
-setattr(vis, "isc", vis_data["isc"])
 
 
 ###############################################################################
@@ -93,6 +95,24 @@
 ###############################################################################
 # MEM
 
+# Compute percent_lambda
+# Loop through ISCs starting with 6-9, but if we don't have at least 2 vis, lower isc_min to include next one down, etc.
+isc_min = 6
+nbig = 0
+
+while isc_min >= 0 and nbig < 2:
+    ibig = np.argwhere(vis.meta.vis_labels >= isc_min)
+    nbig = len(ibig)
+    isc_min = isc_min - 1
+
+# If still don't have at least 2 vis, return -1, otherwise calculate mean (but reject points > sigma away from mean)
+if nbig < 2:
+    snr_value = -1
+else:
+    snr_value, _ = resistant_mean((np.abs(vis.visibilities[ibig]) / vis.amplitude_uncertainty[ibig]).flatten(), 3)
+
+percent_lambda = 11.0 / (snr_value**2 + 383.0)
+
 mem_map = mem(vis, shape=[129, 129] * apu.pixel, pixel_size=[2, 2] * apu.arcsec / apu.pix)
 mem_map.plot()
 

diff --git a/examples/stix.py b/examples/stix.py
@@ -8,26 +8,27 @@
 """
 
 import pickle
-import urllib.request
 
 import astropy.units as apu
 import matplotlib.pyplot as plt
+import numpy as np
 
 from xrayvision.clean import vis_clean
 from xrayvision.imaging import vis_psf_map, vis_to_map
-from xrayvision.mem import mem
+from xrayvision.mem import mem, resistant_mean
 
 ###############################################################################
 # Create images from STIX visibility data.
 #
 # The STIX data has already been prepared and stored in python pickle format
 # the variables can be simply restored.
 
-stix_data = pickle.load(urllib.request.urlopen("https://pub099.cs.technik.fhnw.ch/demo/stix_vis.pkl"))
+with open("./stix_vis.pkl", "rb") as file:
+    stix_data = pickle.load(file)
 
-time_range, energy_range, offset, stix_vis = stix_data
-stix_vis.phase_centre = [0, 0] * apu.arcsec
-stix_vis.offset = offset
+time_range = stix_data["time_range"]
+energy_range = stix_data["energy_range"]
+stix_vis = stix_data["stix_visibilities"]
 
 ###############################################################################
 # Lets have a look at the point spread function (PSF) or dirty beam
@@ -56,7 +57,13 @@
 ###############################################################################
 # MEM
 
-mem_map = mem(stix_vis, shape=[129, 129] * apu.pixel, pixel_size=[2, 2] * apu.arcsec / apu.pix)
+# Compute percent_lambda
+snr_value, _ = resistant_mean((np.abs(stix_vis.visibilities) / stix_vis.amplitude_uncertainty).flatten(), 3)
+percent_lambda = 2 / (snr_value**2 + 90)
+
+mem_map = mem(
+    stix_vis, shape=[129, 129] * apu.pixel, pixel_size=[2, 2] * apu.arcsec / apu.pix, percent_lambda=percent_lambda
+)
 mem_map.plot()
 
 ###############################################################################

diff --git a/examples/stix_vis.pkl b/examples/stix_vis.pkl
diff --git a/xrayvision/clean.py b/xrayvision/clean.py
@@ -21,7 +21,7 @@
 from sunpy.map.map_factory import Map
 
 from xrayvision.imaging import vis_psf_image, vis_to_map
-from xrayvision.visibility import Visibility
+from xrayvision.visibility import Visibilities
 
 __all__ = ["clean", "vis_clean", "ms_clean", "vis_ms_clean"]
 
@@ -98,7 +98,7 @@ def clean(
     #     raise ValueError('')
     pad = [0 if x % 2 == 0 else 1 for x in dirty_map.shape]
 
-    # Assume beam, map phase_centre is in middle
+    # Assume beam, map phase_center is in middle
     beam_center = (dirty_beam.shape[0] - 1) / 2.0, (dirty_beam.shape[1] - 1) / 2.0
     map_center = (dirty_map.shape[0] - 1) / 2.0, (dirty_map.shape[1] - 1) / 2.0
 
@@ -173,7 +173,7 @@ def clean(
 
 @u.quantity_input
 def vis_clean(
-    vis: Visibility,
+    vis: Visibilities,
     shape: Quantity[u.pix],
     pixel_size: Quantity[u.arcsec / u.pix],
     clean_beam_width: Optional[Quantity[u.arcsec]] = 4.0,
@@ -403,7 +403,7 @@ def ms_clean(
 
 
 def vis_ms_clean(
-    vis: Visibility,
+    vis: Visibilities,
     shape: Quantity[u.pix],
     pixel_size: Quantity[u.arcsec / u.pix],
     scales: Optional[Iterable],

diff --git a/xrayvision/imaging.py b/xrayvision/imaging.py
@@ -6,7 +6,7 @@
 from sunpy.map import GenericMap, Map
 
 from xrayvision.transform import dft_map, idft_map
-from xrayvision.visibility import Visibility
+from xrayvision.visibility import Visibilities
 
 __all__ = [
     "get_weights",
@@ -24,7 +24,7 @@
 WEIGHT_SCHEMES = ("natural", "uniform")
 
 
-def get_weights(vis: Visibility, scheme: str = "natural", norm: bool = True) -> np.ndarray:
+def get_weights(vis: Visibilities, scheme: str = "natural", norm: bool = True) -> np.ndarray:
     r"""
     Return spatial frequency weight factors for each visibility.
 
@@ -47,7 +47,7 @@ def get_weights(vis: Visibility, scheme: str = "natural", norm: bool = True) ->
         raise ValueError(f"Invalid weighting scheme {scheme}, must be one of: {WEIGHT_SCHEMES}")
     weights = np.sqrt(vis.u**2 + vis.v**2).value
     if scheme == "natural":
-        weights = np.ones_like(vis.vis, dtype=float)
+        weights = np.ones_like(vis.visibilities, dtype=float)
 
     if norm:
         weights /= weights.sum()
@@ -97,11 +97,11 @@ def image_to_vis(
     *,
     u: Quantity[apu.arcsec**-1],
     v: Quantity[apu.arcsec**-1],
-    phase_centre: Optional[Quantity[apu.arcsec]] = (0.0, 0.0) * apu.arcsec,
+    phase_center: Optional[Quantity[apu.arcsec]] = (0.0, 0.0) * apu.arcsec,
     pixel_size: Optional[Quantity[apu.arcsec / apu.pix]] = 1.0 * apu.arcsec / apu.pix,
-) -> Visibility:
+) -> Visibilities:
     r"""
-    Return a Visibility created from the image and u, v sampling.
+    Return a Visibilities object created from the image and u, v sampling.
 
     Parameters
     ----------
@@ -111,27 +111,29 @@ def image_to_vis(
         Array of u coordinates where the visibilities will be evaluated
     v :
         Array of v coordinates where the visibilities will be evaluated
-    phase_centre :
-        The coordinates the phase_centre.
+    phase_center :
+        The coordinates the phase_center.
     pixel_size :
         Size of pixels, if only one value is passed, assume square pixels (repeating the value).
 
     Returns
     -------
     :
-        The new visibility object
+        The new Visibilities object
 
     """
     pixel_size = validate_and_expand_kwarg(pixel_size, "pixel_size")
     if not (apu.get_physical_type((1 / u).unit) == ANGLE and apu.get_physical_type((1 / v).unit) == ANGLE):
         raise ValueError("u and v must be inverse angle (e.g. 1/deg or 1/arcsec")
-    vis = dft_map(image, u=u, v=v, phase_centre=phase_centre, pixel_size=pixel_size)
-    return Visibility(vis, u=u, v=v, offset=phase_centre)
+    vis = dft_map(
+        image, u=u, v=v, phase_center=[0.0, 0.0] * apu.arcsec, pixel_size=pixel_size
+    )  # TODO: adapt to generic map center
+    return Visibilities(vis, u=u, v=v, phase_center=phase_center)
 
 
 @apu.quantity_input()
 def vis_to_image(
-    vis: Visibility,
+    vis: Visibilities,
     shape: Quantity[apu.pix] = (65, 65) * apu.pixel,
     pixel_size: Optional[Quantity[apu.arcsec / apu.pix]] = 1 * apu.arcsec / apu.pix,
     scheme: str = "natural",
@@ -161,15 +163,21 @@ def vis_to_image(
     shape = shape.to(apu.pixel)
     weights = get_weights(vis, scheme=scheme)
     bp_arr = idft_map(
-        vis.vis, u=vis.u, v=vis.v, shape=shape, weights=weights, pixel_size=pixel_size, phase_centre=vis.phase_centre
+        vis.visibilities,
+        u=vis.u,
+        v=vis.v,
+        shape=shape,
+        weights=weights,
+        pixel_size=pixel_size,
+        phase_center=[0.0, 0.0] * apu.arcsec,  # TODO update to have generic image center
     )
 
     return bp_arr
 
 
 @apu.quantity_input
 def vis_psf_map(
-    vis: Visibility,
+    vis: Visibilities,
     *,
     shape: Quantity[apu.pix] = (65, 65) * apu.pixel,
     pixel_size: Optional[Quantity[apu.arcsec / apu.pix]] = 1 * apu.arcsec / apu.pix,
@@ -203,7 +211,7 @@ def vis_psf_map(
 
 @apu.quantity_input()
 def vis_psf_image(
-    vis: Visibility,
+    vis: Visibilities,
     *,
     shape: Quantity[apu.pix] = (65, 65) * apu.pixel,
     pixel_size: Quantity[apu.arcsec / apu.pix] = 1 * apu.arcsec / apu.pix,
@@ -237,14 +245,19 @@ def vis_psf_image(
     # Make sure psf is always odd so power is in exactly one pixel
     shape = [s // 2 * 2 + 1 for s in shape.to_value(apu.pix)] * shape.unit
     psf_arr = idft_map(
-        np.ones(vis.vis.shape) * vis.vis.unit, u=vis.u, v=vis.v, shape=shape, weights=weights, pixel_size=pixel_size
+        np.ones(vis.visibilities.shape) * vis.visibilities.unit,
+        u=vis.u,
+        v=vis.v,
+        shape=shape,
+        weights=weights,
+        pixel_size=pixel_size,
     )
     return psf_arr
 
 
 @apu.quantity_input()
 def vis_to_map(
-    vis: Visibility,
+    vis: Visibilities,
     shape: Quantity[apu.pix] = (65, 65) * apu.pixel,
     pixel_size: Optional[Quantity[apu.arcsec / apu.pix]] = 1 * apu.arcsec / apu.pixel,
     scheme: Optional[str] = "natural",
@@ -278,7 +291,7 @@ def vis_to_map(
 
 
 @apu.quantity_input()
-def generate_header(vis: Visibility, *, shape: Quantity[apu.pix], pixel_size: Quantity[apu.arcsec / apu.pix]) -> dict:
+def generate_header(vis: Visibilities, *, shape: Quantity[apu.pix], pixel_size: Quantity[apu.arcsec / apu.pix]) -> dict:
     r"""
     Generate a map head given the visibilities, pixel size and shape
 
@@ -296,8 +309,8 @@ def generate_header(vis: Visibility, *, shape: Quantity[apu.pix], pixel_size: Qu
     :
     """
     header = {
-        "crval1": (vis.offset[1]).to_value(apu.arcsec),
-        "crval2": (vis.offset[0]).to_value(apu.arcsec),
+        "crval1": (vis.phase_center[1]).to_value(apu.arcsec),
+        "crval2": (vis.phase_center[0]).to_value(apu.arcsec),
         "cdelt1": (pixel_size[1] * apu.pix).to_value(apu.arcsec),
         "cdelt2": (pixel_size[0] * apu.pix).to_value(apu.arcsec),
         "ctype1": "HPLN-TAN",
@@ -312,9 +325,9 @@ def generate_header(vis: Visibility, *, shape: Quantity[apu.pix], pixel_size: Qu
 
 
 @apu.quantity_input()
-def map_to_vis(amap: GenericMap, *, u: Quantity[1 / apu.arcsec], v: Quantity[1 / apu.arcsec]) -> Visibility:
+def map_to_vis(amap: GenericMap, *, u: Quantity[1 / apu.arcsec], v: Quantity[1 / apu.arcsec]) -> Visibilities:
     r"""
-    Return a Visibility object created from the map, sampling it at give `u`, `v` coordinates.
+    Return a Visibilities object created from the map, sampling it at give `u`, `v` coordinates.
 
     Parameters
     ----------
@@ -328,7 +341,7 @@ def map_to_vis(amap: GenericMap, *, u: Quantity[1 / apu.arcsec], v: Quantity[1 /
     Returns
     -------
     :
-        The new visibility object
+        The new Visibilities object
 
     """
     if not apu.get_physical_type(1 / u) == ANGLE and apu.get_physical_type(1 / v) == ANGLE:
@@ -346,6 +359,8 @@ def map_to_vis(amap: GenericMap, *, u: Quantity[1 / apu.arcsec], v: Quantity[1 /
     if "cdelt2" in meta:
         new_psize[0] = float(meta["cdelt2"])
 
-    vis = image_to_vis(amap.data, u=u, v=v, pixel_size=new_psize * apu.arcsec / apu.pix)
-    vis.offset = new_pos * apu.arcsec
+    vis = image_to_vis(
+        amap.quantity, u=u, v=v, pixel_size=new_psize * apu.arcsec / apu.pix, phase_center=new_pos * apu.arcsec
+    )
+
     return vis
Original file line number	Diff line number	Diff line change
		@@ -0,0 +1 @@
		Make the software compatible with :class:`~xrayvision.visibility.Visibilities`.