Update for new Visibilities object and some tidy

changelog/124.feature.rst

@@ -0,0 +1 @@
+Update code and examples to use new :class:`~xrayvision.visibility.Visibilities`, tidy API of `~stixpy.imaging.em.em`.

docs/reference/index.rst

@@ -11,6 +11,7 @@ Reference
    coordinates
    data
    map
+   net
    product
    science
    timeseries

docs/reference/net.rst

@@ -0,0 +1,12 @@
+.. _net:
+
+Net ('stixpy.net')
+******************
+
+The ``net`` submodule contains a Fido client and attrs for STIX data
+
+.. automodapi:: stixpy.net
+
+.. automodapi:: stixpy.net.attrs
+
+.. automodapi:: stixpy.net.client

examples/imaging_demo.py

@@ -16,20 +16,14 @@
 import astropy.units as u
 import matplotlib.pyplot as plt
 import numpy as np
-from astropy.time import Time
 from sunpy.coordinates import HeliographicStonyhurst, Helioprojective
 from sunpy.map import Map, make_fitswcs_header
+from sunpy.time import TimeRange
 from xrayvision.clean import vis_clean
 from xrayvision.imaging import vis_to_image, vis_to_map
 from xrayvision.mem import mem, resistant_mean
-from xrayvision.visibility import Visibilities, VisMeta
 
-from stixpy.calibration.visibility import (
-    calibrate_visibility,
-    create_meta_pixels,
-    create_visibility,
-    get_uv_points_data,
-)
+from stixpy.calibration.visibility import calibrate_visibility, create_meta_pixels, create_visibility
 from stixpy.coordinates.frames import STIXImaging
 from stixpy.coordinates.transforms import get_hpc_info
 from stixpy.imaging.em import em
@@ -62,57 +56,51 @@
     "2021-09-23T09:00:00",
     "2021-09-23T12:00:00",
 ]  # Set this range larger than the actual observation time
-energy_range = [28, 40]
+energy_range = [28, 40] * u.keV
 
 ###############################################################################
 # Create the meta pixel, A, B, C, D for the science and the background data
 
 meta_pixels_sci = create_meta_pixels(
-    cpd_sci, time_range=time_range_sci, energy_range=energy_range, phase_center=[0, 0] * u.arcsec, no_shadowing=True
+    cpd_sci, time_range=time_range_sci, energy_range=energy_range, flare_location=[0, 0] * u.arcsec, no_shadowing=True
 )
 
 meta_pixels_bkg = create_meta_pixels(
-    cpd_bkg, time_range=time_range_bkg, energy_range=energy_range, phase_center=[0, 0] * u.arcsec, no_shadowing=True
+    cpd_bkg, time_range=time_range_bkg, energy_range=energy_range, flare_location=[0, 0] * u.arcsec, no_shadowing=True
 )
 
 ###############################################################################
 # Perform background subtraction
 
 meta_pixels_bkg_subtracted = {
+    **meta_pixels_sci,
     "abcd_rate_kev_cm": meta_pixels_sci["abcd_rate_kev_cm"] - meta_pixels_bkg["abcd_rate_kev_cm"],
     "abcd_rate_error_kev_cm": np.sqrt(
         meta_pixels_sci["abcd_rate_error_kev_cm"] ** 2 + meta_pixels_bkg["abcd_rate_error_kev_cm"] ** 2
     ),
 }
 
 ###############################################################################
-# Create visibilites from the meta pixels
+# Create visibilities from the meta pixels
 
 vis = create_visibility(meta_pixels_bkg_subtracted)
 
 ###############################################################################
 # Calibrate the visibilities
 
-# Extra phase calibration not needed with these
-uv_data = get_uv_points_data()
-vis.u = uv_data["u"]
-vis.v = uv_data["v"]
-
 cal_vis = calibrate_visibility(vis)
 
 ###############################################################################
 # Selected detectors 10 to 7
 
 # order by sub-collimator e.g 10a, 10b, 10c, 9a, 9b, 9c ....
 isc_10_7 = [3, 20, 22, 16, 14, 32, 21, 26, 4, 24, 8, 28]
-idx = np.argwhere(np.isin(cal_vis.isc, isc_10_7)).ravel()
+idx = np.argwhere(np.isin(cal_vis.meta["isc"], isc_10_7)).ravel()
 
 ###############################################################################
-# Create an ``xrayvsion`` Visibilities object
+# Slice the visibilities to detectors 10 - 7
 
-stix_vis = Visibilities(
-    cal_vis.obsvis[idx], cal_vis.u[idx], cal_vis.v[idx], amplitude_uncertainty=cal_vis.amplitude_error[idx]
-)
+vis10_7 = cal_vis[idx]
 
 ###############################################################################
 # Set up image parameters
@@ -123,23 +111,39 @@
 ###############################################################################
 # Make a full disk back projection (inverse transform) map
 
-bp_image = vis_to_image(stix_vis, imsize, pixel_size=pixel)
+bp_image = vis_to_image(vis10_7, imsize, pixel_size=pixel)
 
-date_avg = Time("2021-09-23T15:22:30")
-roll, solo_xyz, pointing = get_hpc_info(date_avg)
+###############################################################################
+# Obtain the necessary metadata to create a sunpy map in the STIXImaging frame
 
-solo = HeliographicStonyhurst(*solo_xyz, obstime=date_avg, representation_type="cartesian")
-coord = STIXImaging(0 * u.arcsec, 0 * u.arcsec, obstime="2021-09-23T15:22:30", observer=solo)
+vis_tr = TimeRange(vis.meta["time_range"])
+roll, solo_xyz, pointing = get_hpc_info(vis_tr.start, vis_tr.end)
+solo = HeliographicStonyhurst(*solo_xyz, obstime=vis_tr.center, representation_type="cartesian")
+coord = STIXImaging(0 * u.arcsec, 0 * u.arcsec, obstime=vis_tr.start, observer=solo)
 header = make_fitswcs_header(
     bp_image, coord, telescope="STIX", observatory="Solar Orbiter", scale=[10, 10] * u.arcsec / u.pix
 )
 fd_bp_map = Map((bp_image, header))
 
+###############################################################################
+# Convert the coordinates and make a map in Helioprojective and rotate so "North" is "up"
+
 hpc_ref = coord.transform_to(Helioprojective(observer=solo, obstime=fd_bp_map.date))  # Center of STIX pointing in HPC
 header_hp = make_fitswcs_header(bp_image, hpc_ref, scale=[10, 10] * u.arcsec / u.pix, rotation_angle=90 * u.deg + roll)
 hp_map = Map((bp_image, header_hp))
 hp_map_rotated = hp_map.rotate()
 
+###############################################################################
+# Estimate the flare location and plot on top of back projection map. Note the coordinates
+# are automaiccally converted from the STIXImaging to Helioprojective
+
+max_pixel = np.argwhere(fd_bp_map.data == fd_bp_map.data.max()).ravel() * u.pixel
+# because WCS axes and array are reversed
+max_stix = fd_bp_map.pixel_to_world(max_pixel[1], max_pixel[0])
+
+###############################################################################
+# Plot the both maps
+
 fig = plt.figure(figsize=(12, 8))
 ax0 = fig.add_subplot(1, 2, 1, projection=fd_bp_map)
 ax1 = fig.add_subplot(1, 2, 2, projection=hp_map_rotated)
@@ -151,58 +155,38 @@
 hp_map_rotated.draw_limb()
 hp_map_rotated.draw_grid()
 
-###############################################################################
-# Estimate the flare location and plot on top of back projection map.
-
-max_pixel = np.argwhere(fd_bp_map.data == fd_bp_map.data.max()).ravel() * u.pixel
-# because WCS axes and array are reversed
-max_stix = fd_bp_map.pixel_to_world(max_pixel[1], max_pixel[0])
-
 ax0.plot_coord(max_stix, marker=".", markersize=50, fillstyle="none", color="r", markeredgewidth=2)
 ax1.plot_coord(max_stix, marker=".", markersize=50, fillstyle="none", color="r", markeredgewidth=2)
 
-
 ################################################################################
 # Use estimated flare location to create more accurate visibilities
 
-# because WCS axes and array are reversed and all positions are expected follow array indices
-flare_pos_stix = [max_stix.Ty.to_value("arcsec"), max_stix.Tx.to_value("arcsec")] * u.arcsec
-
 meta_pixels_sci = create_meta_pixels(
-    cpd_sci, time_range=time_range_sci, energy_range=energy_range, phase_center=flare_pos_stix, no_shadowing=True
+    cpd_sci, time_range=time_range_sci, energy_range=energy_range, flare_location=max_stix, no_shadowing=True
 )
 
 meta_pixels_bkg_subtracted = {
+    **meta_pixels_sci,
     "abcd_rate_kev_cm": meta_pixels_sci["abcd_rate_kev_cm"] - meta_pixels_bkg["abcd_rate_kev_cm"],
     "abcd_rate_error_kev_cm": np.sqrt(
         meta_pixels_sci["abcd_rate_error_kev_cm"] ** 2 + meta_pixels_bkg["abcd_rate_error_kev_cm"] ** 2
     ),
 }
 
 vis = create_visibility(meta_pixels_bkg_subtracted)
-uv_data = get_uv_points_data()
-vis.u = uv_data["u"]
-vis.v = uv_data["v"]
-cal_vis = calibrate_visibility(vis, flare_pos_stix)
+cal_vis = calibrate_visibility(vis, flare_location=max_stix)
 
 ###############################################################################
 # Selected detectors 10 to 3
 # order by sub-collimator e.g 10a, 10b, 10c, 9a, 9b, 9c ....
 isc_10_3 = [3, 20, 22, 16, 14, 32, 21, 26, 4, 24, 8, 28, 15, 27, 31, 6, 30, 2, 25, 5, 23, 7, 29, 1]
-idx = np.argwhere(np.isin(cal_vis.isc, isc_10_3)).ravel()
+idx = np.argwhere(np.isin(cal_vis.meta["isc"], isc_10_3)).ravel()
 
 ###############################################################################
 # Create an ``xrayvsion`` visibility object
 
-meta = VisMeta({"vis_labels": cal_vis.isc[idx]})
-stix_vis1 = Visibilities(
-    cal_vis.obsvis[idx],
-    cal_vis.u[idx],
-    cal_vis.v[idx],
-    amplitude_uncertainty=cal_vis.amplitude_error[idx],
-    phase_center=flare_pos_stix,
-    meta=meta,
-)
+cal_vis.meta["offset"] = max_stix
+vis10_3 = cal_vis[idx]
 
 ###############################################################################
 # Set up image parameters
@@ -213,52 +197,68 @@
 ###############################################################################
 # Create a back projection image with natural weighting
 
-bp_nat = vis_to_image(stix_vis1, imsize, pixel_size=pixel)
+bp_nat = vis_to_image(vis10_3, imsize, pixel_size=pixel)
 
 ###############################################################################
 # Create a back projection image with uniform weighting
 
-bp_uni = vis_to_image(stix_vis1, imsize, pixel_size=pixel, scheme="uniform")
+bp_uni = vis_to_image(vis10_3, imsize, pixel_size=pixel, scheme="uniform")
 
 ###############################################################################
 # Create a `sunpy.map.Map` with back projection
 
-bp_map = vis_to_map(stix_vis1, imsize, pixel_size=pixel)
+bp_map = vis_to_map(vis10_3, imsize, pixel_size=pixel)
 
 ###############################################################################
-# Crete a map using the clean algorithm `vis_clean`
+# Crete a clean image using the clean algorithm `vis_clean`
 
 niter = 200  # number of iterations
 gain = 0.1  # gain used in each clean iteration
 beam_width = 20.0 * u.arcsec
 clean_map, model_map, resid_map = vis_clean(
-    stix_vis1, imsize, pixel_size=pixel, gain=gain, niter=niter, clean_beam_width=20 * u.arcsec
+    vis10_3, imsize, pixel_size=pixel, gain=gain, niter=niter, clean_beam_width=20 * u.arcsec
 )
 
+###############################################################################
+# Create a sunpy map for the clean image in Helioprojective
+
+header = make_fitswcs_header(
+    clean_map.data,
+    max_stix.transform_to(Helioprojective(obstime=vis_tr.center, observer=solo)),
+    telescope="STIX",
+    observatory="Solar Orbiter",
+    scale=pixel,
+    rotation_angle=90 * u.deg + roll,
+)
+
+clean_map = Map((clean_map.data, header))
+plt.figure()
+clean_map.rotate().plot()
+
 ###############################################################################
 # Crete a map using the MEM GE algorithm `mem`
 
-snr_value, _ = resistant_mean((np.abs(stix_vis1.visibilities) / stix_vis1.amplitude_uncertainty).flatten(), 3)
+snr_value, _ = resistant_mean((np.abs(vis10_3.visibilities) / vis10_3.amplitude_uncertainty).flatten(), 3)
 percent_lambda = 2 / (snr_value**2 + 90)
-mem_map = mem(stix_vis1, shape=imsize, pixel_size=pixel, percent_lambda=percent_lambda)
+mem_map = mem(vis10_3, shape=imsize, pixel_size=pixel, percent_lambda=percent_lambda)
 
 
 ###############################################################################
 # Crete a map using the EM algorithm `EM`
 
 em_map = em(
     meta_pixels_bkg_subtracted["abcd_rate_kev_cm"],
-    stix_vis1,
+    cal_vis,
     shape=imsize,
     pixel_size=pixel,
-    flare_xy=flare_pos_stix,
+    flare_location=max_stix,
     idx=idx,
 )
 
 vmax = max([clean_map.data.max(), mem_map.data.max(), em_map.value.max()])
 
 ###############################################################################
-# Compare the image from each algorithm
+# Finally compare the images from each algorithm
 
 fig, axes = plt.subplots(2, 2)
 a = axes[0, 0].imshow(bp_nat.value)
@@ -275,9 +275,3 @@
 fig.colorbar(d)
 fig.tight_layout()
 plt.show()
-
-# roll, pos, ref = get_hpc_info(Time(time_range[0]))
-# solo_coord = SkyCoord(*pos, frame='heliocentricearthecliptic', representation_type='cartesian', obstime=Time(time_range[0]))
-# ref_coord = SkyCoord(ref[0], ref[1], obstime=Time(time_range[0]), observer=solo_coord, frame='helioprojective')
-# header = make_fitswcs_header(fd_bp_map.value, ref_coord, rotation_angle=-roll+90*u.deg, scale=[10,10]*u.arcsec/u.pix)
-# stix_hpc_map = Map(fd_bp_map.value, header)

setup.cfg

@@ -36,7 +36,7 @@ install_requires =
     intervaltree
     roentgen
     matplotlib
-    xrayvisim @ git+https://github.com/TCDSolar/xrayvision.git
+    xrayvisim~=0.2.0rc
 [options.extras_require]
 test =
     pytest

stixpy/calibration/tests/test_visibility.py

@@ -1,8 +1,10 @@
 import astropy.units as u
 import pytest
 from astropy.tests.helper import assert_quantity_allclose
+from astropy.time import Time
 
 from stixpy.calibration.visibility import create_meta_pixels, get_uv_points_data
+from stixpy.coordinates.frames import STIXImaging
 from stixpy.product import Product
 
 
@@ -33,13 +35,13 @@ def test_get_uv_points_data():
 
 
 def test_create_meta_pixels(background_cpd):
-    time_range = ["2022-08-24T14:00:37.271", "2022-08-24T14:50:17.271"]
-    energy_range = [20, 76]
+    time_range = Time(["2022-08-24T14:00:37.271", "2022-08-24T14:50:17.271"])
+    energy_range = [20, 76] * u.keV
     meta_pixels = create_meta_pixels(
         background_cpd,
         time_range=time_range,
         energy_range=energy_range,
-        phase_center=[0, 0] * u.arcsec,
+        flare_location=STIXImaging(0 * u.arcsec, 0 * u.arcsec),
         no_shadowing=True,
     )