From 4d7b12cf028755ad628a4768c092b826864ca9b6 Mon Sep 17 00:00:00 2001 From: Shane Maloney Date: Wed, 15 May 2024 11:36:44 +0100 Subject: [PATCH] Updates for latest xrayvision --- examples/imaging_demo.py | 29 ++++++++++++++++------------- stixpy/imaging/em.py | 9 +++------ 2 files changed, 19 insertions(+), 19 deletions(-) diff --git a/examples/imaging_demo.py b/examples/imaging_demo.py index c5b0ed1..bf8e545 100644 --- a/examples/imaging_demo.py +++ b/examples/imaging_demo.py @@ -114,7 +114,7 @@ # Set up image parameters imsize = [512, 512] * u.pixel # number of pixels of the map to reconstruct -pixel = [10, 10] * u.arcsec # pixel size in aresec +pixel = [10, 10] * u.arcsec/u.pixel # pixel size in aresec ############################################################################### # Make a full disk back projection (inverse transform) map @@ -149,18 +149,21 @@ # 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 are reverse order -max_hpc = fd_bp_map.pixel_to_world(max_pixel[1], max_pixel[0]) +# 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_hpc, marker=".", markersize=50, fillstyle="none", color="r", markeredgewidth=2) -ax1.plot_coord(max_hpc, marker=".", markersize=50, fillstyle="none", color="r", markeredgewidth=2) +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=[max_hpc.Tx, max_hpc.Ty], no_shadowing=True + cpd_sci, time_range=time_range_sci, energy_range=energy_range, phase_center=flare_pos_stix, no_shadowing=True ) meta_pixels_bkg_subtracted = {"abcd_rate_kev_cm": meta_pixels_sci["abcd_rate_kev_cm"] - meta_pixels_bkg["abcd_rate_kev_cm"], @@ -171,7 +174,7 @@ uv_data = get_uv_points_data() vis.u = uv_data['u'] vis.v = uv_data['v'] -cal_vis = calibrate_visibility(vis, [max_hpc.Tx, max_hpc.Ty]) +cal_vis = calibrate_visibility(vis, flare_pos_stix) ############################################################################### # Selected detectors 10 to 3 @@ -186,7 +189,7 @@ vis=cal_vis.obsvis[idx], u=cal_vis.u[idx], v=cal_vis.v[idx], - offset=np.array([max_hpc.Tx.value, max_hpc.Ty.value]) * u.arcsec, + offset=flare_pos_stix, ) skeys = stix_vis1.__dict__.keys() for k, v in cal_vis.__dict__.items(): @@ -197,7 +200,7 @@ # Set up image parameters imsize = [129, 129] * u.pixel # number of pixels of the map to reconstruct -pixel = [2, 2] * u.arcsec # pixel size in arcsec +pixel = [2, 2] * u.arcsec / u.pixel # pixel size in arcsec ############################################################################### # Create a back projection image with natural weighting @@ -207,7 +210,7 @@ ############################################################################### # Create a back projection image with uniform weighting -bp_uni = vis_to_image(stix_vis1, imsize, pixel_size=pixel, natural=False) +bp_uni = vis_to_image(stix_vis1, imsize, pixel_size=pixel, scheme='uniform') ############################################################################### # Create a `sunpy.map.Map` with back projection @@ -221,7 +224,7 @@ 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=pixel, gain=gain, niter=niter, clean_beam_width=20 * u.arcsec + stix_vis1, imsize, pixel_size=pixel, gain=gain, niter=niter, clean_beam_width=20 * u.arcsec ) ############################################################################### @@ -237,7 +240,7 @@ vis=cal_vis.obsvis[idx], u=cal_vis.u[idx], v=cal_vis.v[idx], - center=np.array([max_hpc.Tx.value, max_hpc.Ty.value]) * u.arcsec, + offset=flare_pos_stix, ) skeys = stix_vis1.__dict__.keys() for k, v in cal_vis.__dict__.items(): @@ -249,7 +252,7 @@ stix_vis1, shape=imsize, pixel_size=pixel, - flare_xy=np.array([max_hpc.Tx.value, max_hpc.Ty.value]) * u.arcsec, + flare_xy=flare_pos_stix, idx=idx, ) diff --git a/stixpy/imaging/em.py b/stixpy/imaging/em.py index 6f5742e..97e4a1f 100644 --- a/stixpy/imaging/em.py +++ b/stixpy/imaging/em.py @@ -58,12 +58,9 @@ def get_transmission_matrix(u, v, shape=[64, 64] * apu.pix, pixel_size=[4.0, 4.0 idx = [6, 28, 0, 24, 4, 22, 5, 29, 1, 14, 26, 30, 23, 7, 27, 20, 25, 3, 15, 13, 31, 2, 19, 21] phase_cor = phase_cor[idx] - m, n = shape.to_value("pix") + x = generate_xy(shape[0], pixel_size=pixel_size[0], phase_centre=center[0]) + y = generate_xy(shape[1], pixel_size=pixel_size[1], phase_centre=center[1]) - # center = center - [26.1, 58.2] * apu.arcsec - - y = generate_xy(m, pixel_size=pixel_size[1], center=center[1]) - x = generate_xy(n, pixel_size=pixel_size[0], center=center[0]) x, y = np.meshgrid(x, y) # Check apu are correct for exp need to be dimensionless and then remove apu for speed if (u[0] * x[0, 0]).unit == apu.dimensionless_unscaled and (v[0] * y[0, 0]).unit == apu.dimensionless_unscaled: @@ -154,5 +151,5 @@ def em(countrates, vis, shape, pixel_size, maxiter=5000, tolerance=0.001, *, fla m, n = shape.to_value("pix").astype(int) # Note the order keyword here to match `y.flatten(order='F')` above at line 108 # TODO investigate why .T is needed probably because of fortran ordering - x_im = x.reshape(m, n, order="F").T + x_im = x.reshape(m, n, order="F") return x_im