From 540cf707846cb698f061b60e8916bc587108d05d Mon Sep 17 00:00:00 2001 From: Victor de Souza Magalhaes Date: Tue, 18 Feb 2025 12:03:43 -0700 Subject: [PATCH 01/53] AntennaSurface class now passes the dish outer limit rather than its maximal radius to create_aperture mask. --- src/astrohack/antenna/antenna_surface.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/astrohack/antenna/antenna_surface.py b/src/astrohack/antenna/antenna_surface.py index 081c02a4..f1aac746 100644 --- a/src/astrohack/antenna/antenna_surface.py +++ b/src/astrohack/antenna/antenna_surface.py @@ -303,7 +303,7 @@ def _create_aperture_mask(self, clip_type, clip_level, exclude_shadows): arm_angle = 0.0 self.base_mask, self.rad, self.phi = create_aperture_mask(self.u_axis, self.v_axis, self.telescope.inlim, - self.telescope.diam/2, + self.telescope.oulim, arm_width=arm_width, arm_angle=arm_angle, return_polar_meshes=True) From 4d60cf2f8c0ae4f6c6b975ec26da4d81d802a504 Mon Sep 17 00:00:00 2001 From: Victor de Souza Magalhaes Date: Tue, 18 Feb 2025 12:04:47 -0700 Subject: [PATCH 02/53] Implemented support for arm shadows with different widths for different radii ranges --- src/astrohack/utils/algorithms.py | 38 +++++++++++++++++++++---------- 1 file changed, 26 insertions(+), 12 deletions(-) diff --git a/src/astrohack/utils/algorithms.py b/src/astrohack/utils/algorithms.py index a6895e9a..5e52193c 100644 --- a/src/astrohack/utils/algorithms.py +++ b/src/astrohack/utils/algorithms.py @@ -547,21 +547,35 @@ def create_aperture_mask(x_axis, y_axis, inner_rad, outer_rad, arm_width=None, a if arm_width is None: pass + elif isinstance(arm_width, (float, int)): + mask = _arm_shadow_masking(mask, x_mesh, y_mesh, radius_mesh, inner_rad, outer_rad, arm_width, arm_angle) + elif isinstance(arm_width, list): + for section in arm_width: + minradius, maxradius, width = section + mask = _arm_shadow_masking(mask, x_mesh, y_mesh, radius_mesh, minradius, maxradius, width, arm_angle) + else: - if arm_angle % pi/2 == 0: - mask = np.where(np.abs(x_mesh) < arm_width/2., False, mask) - mask = np.where(np.abs(y_mesh) < arm_width/2., False, mask) - else: - # first shadow - coeff = np.tan(arm_angle % pi) - distance = np.abs((coeff*x_mesh-y_mesh)/np.sqrt(coeff**2+1)) - mask = np.where(distance < arm_width/2., False, mask) - # second shadow - coeff = np.tan(arm_angle % pi + pi/2) - distance = np.abs((coeff*x_mesh-y_mesh)/np.sqrt(coeff**2+1)) - mask = np.where(distance < arm_width/2., False, mask) + raise Exception(f"Don't know how to handle an arm width of class {type(arm_width)}") if return_polar_meshes: return mask, radius_mesh, polar_angle_mesh else: return mask + + +def _arm_shadow_masking(inmask, x_mesh, y_mesh, radius_mesh, minradius, maxradius, width, angle): + radial_mask = np.where(radius_mesh < minradius, False, inmask) + radial_mask = np.where(radius_mesh >= maxradius, False, radial_mask) + if angle % pi/2 == 0: + oumask = np.where(np.bitwise_and(np.abs(x_mesh) < width/2., radial_mask), False, inmask) + oumask = np.where(np.bitwise_and(np.abs(y_mesh) < width/2., radial_mask), False, oumask) + else: + # first shadow + coeff = np.tan(angle % pi) + distance = np.abs((coeff*x_mesh-y_mesh)/np.sqrt(coeff**2+1)) + oumask = np.where(np.bitwise_and(distance < width/2., radial_mask), False, inmask) + # second shadow + coeff = np.tan(angle % pi + pi/2) + distance = np.abs((coeff*x_mesh-y_mesh)/np.sqrt(coeff**2+1)) + oumask = np.where(np.bitwise_and(distance < width/2., radial_mask), False, oumask) + return oumask From d3134cd969f3c50c3f52bfd27a873c4da1d3807e Mon Sep 17 00:00:00 2001 From: Victor de Souza Magalhaes Date: Tue, 18 Feb 2025 12:08:26 -0700 Subject: [PATCH 03/53] Added minimal documentation to create_aperture_mask in algorithms.py --- src/astrohack/utils/algorithms.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/astrohack/utils/algorithms.py b/src/astrohack/utils/algorithms.py index 5e52193c..bf9e2f47 100644 --- a/src/astrohack/utils/algorithms.py +++ b/src/astrohack/utils/algorithms.py @@ -538,6 +538,18 @@ def create_coordinate_images(x_axis, y_axis, create_polar_coordinates=False): def create_aperture_mask(x_axis, y_axis, inner_rad, outer_rad, arm_width=None, arm_angle=0, return_polar_meshes=False): """ + Create a basic aperture mask with support for feed supporting arms shadows + Args: + x_axis: The X axis of the Aperture + y_axis: The Y axis of the Aperture + inner_rad: The innermost radius for valid data in aperture + outer_rad: The outermost radius for valid data in aperture + arm_width: The width of the feed arm shadows, can be a list with limiting radii or a single value. + arm_angle: The angle between the arm shadows and the X axis + return_polar_meshes: Return the radial and polar meshes to avoid duplicate computations. + + Returns: + """ x_mesh, y_mesh, radius_mesh, polar_angle_mesh = \ create_coordinate_images(x_axis, y_axis, create_polar_coordinates=True) From 3c28c724a3c2fa8c7ec518e69a5642658cd50e2c Mon Sep 17 00:00:00 2001 From: Victor de Souza Magalhaes Date: Tue, 18 Feb 2025 14:19:21 -0700 Subject: [PATCH 04/53] Added a _save_to_dist method to the Telescope class so that modifications made are stored to the astrohack distribution without the problematic keys that cause tests to fail when reading recently modified telescope files. --- src/astrohack/antenna/telescope.py | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/src/astrohack/antenna/telescope.py b/src/astrohack/antenna/telescope.py index 9e2c2dff..79e526dc 100644 --- a/src/astrohack/antenna/telescope.py +++ b/src/astrohack/antenna/telescope.py @@ -110,9 +110,19 @@ def write(self, filename): Args: filename: Name of the output file """ - ledict = vars(self) + obj_dict = vars(self) xds = xr.Dataset() - xds.attrs = ledict + xds.attrs = obj_dict + xds.to_zarr(filename, mode="w", compute=True, consolidated=True) + return + + def _save_to_dist(self): + obj_dict = vars(self) + filename = f'{self.filepath}/{self.filename}' + obj_dict.pop('filepath', None) + obj_dict.pop('filename', None) + xds = xr.Dataset() + xds.attrs = obj_dict xds.to_zarr(filename, mode="w", compute=True, consolidated=True) return @@ -139,7 +149,7 @@ def print(self): def __repr__(self): outstr = '' - ledict = vars(self) - for key, item in ledict.items(): + obj_dict = vars(self) + for key, item in obj_dict.items(): outstr += f"{key:20s} = {str(item)}\n" return outstr From c24a158d9a7018970463ba420807ba74c45a4c16 Mon Sep 17 00:00:00 2001 From: Victor de Souza Magalhaes Date: Tue, 18 Feb 2025 14:19:49 -0700 Subject: [PATCH 05/53] Added staggered arm shadows for the VLA telescope file. --- .../data/telescopes/vla.zarr/.zattrs | 25 +++++++++++++++++-- .../data/telescopes/vla.zarr/.zmetadata | 25 +++++++++++++++++-- 2 files changed, 46 insertions(+), 4 deletions(-) diff --git a/src/astrohack/data/telescopes/vla.zarr/.zattrs b/src/astrohack/data/telescopes/vla.zarr/.zattrs index f063e385..f5dae34c 100644 --- a/src/astrohack/data/telescopes/vla.zarr/.zattrs +++ b/src/astrohack/data/telescopes/vla.zarr/.zattrs @@ -29,8 +29,29 @@ "ea27", "ea28" ], - "arm_shadow_rotation": 0.0, - "arm_shadow_width": 1.5, + "arm_shadow_rotation": 0, + "arm_shadow_width": [ + [ + 1.983, + 7.391, + 0.8 + ], + [ + 7.391, + 9.144, + 1.2 + ], + [ + 9.144, + 10.87, + 1.4 + ], + [ + 10.87, + 12.5, + 1.8 + ] + ], "array_center": { "m0": { "unit": "rad", diff --git a/src/astrohack/data/telescopes/vla.zarr/.zmetadata b/src/astrohack/data/telescopes/vla.zarr/.zmetadata index dbc7a9ac..eae7956c 100644 --- a/src/astrohack/data/telescopes/vla.zarr/.zmetadata +++ b/src/astrohack/data/telescopes/vla.zarr/.zmetadata @@ -31,8 +31,29 @@ "ea27", "ea28" ], - "arm_shadow_rotation": 0.0, - "arm_shadow_width": 1.5, + "arm_shadow_rotation": 0, + "arm_shadow_width": [ + [ + 1.983, + 7.391, + 0.8 + ], + [ + 7.391, + 9.144, + 1.2 + ], + [ + 9.144, + 10.87, + 1.4 + ], + [ + 10.87, + 12.5, + 1.8 + ] + ], "array_center": { "m0": { "unit": "rad", From 1ae4d1ea526db7530e22f238048403f92a91c990 Mon Sep 17 00:00:00 2001 From: Victor de Souza Magalhaes Date: Tue, 18 Feb 2025 14:37:29 -0700 Subject: [PATCH 06/53] Brought in fits comparison CL tool into astrohack, it needs heavy revision before getting a user interface. --- src/astrohack/core/image_compare_tool.py | 362 +++++++++++++++++++++++ 1 file changed, 362 insertions(+) create mode 100644 src/astrohack/core/image_compare_tool.py diff --git a/src/astrohack/core/image_compare_tool.py b/src/astrohack/core/image_compare_tool.py new file mode 100644 index 00000000..17ff0c0e --- /dev/null +++ b/src/astrohack/core/image_compare_tool.py @@ -0,0 +1,362 @@ +from astropy.io import fits +import numpy as np +from scipy.interpolate import griddata +from matplotlib import pyplot as plt +from mpl_toolkits.axes_grid1 import make_axes_locatable +import argparse +from astrohack.visualization.plot_tools import well_positioned_colorbar +from astrohack.visualization.plot_tools import close_figure, get_proper_color_map +import datetime + +parser = argparse.ArgumentParser(description="Compare aperture FTIS maps produced by AIPS and AstroHACK\n" + "Beam maps are not supported") +parser.add_argument('first', type=str, help='first aperture FITS (usually AIPS)') +parser.add_argument('second', type=str, + help='Second aperture FITS (usually AstroHACK)') +parser.add_argument('-n', '--noise-map', action='store_true', + default=False, help='Save noise images') +parser.add_argument('-c', '--noise-clip', type=float, default=1, + help='Noise clipping level (in sigmas)') +parser.add_argument('-q', '--quiet', action='store_true', default=False, + help='Do not print value to screen') +parser.add_argument('-d', '--diameter', type=float, default=25.0, + help='Dish diameter') +parser.add_argument('-b', '--blocage', type=float, default=2.0, + help='Dish inner blocage radius') +parser.add_argument('-m', '--colormap', type=str, default='viridis', + help='Colormap for non residual maps') +parser.add_argument('-p', '--no-division', action='store_true', default=False, + help='Do not perform division, i.e. factor assumed to be 1') +parser.add_argument('-s', '--shadow-width', type=float, default=1.5, + help='Arm shadow width in meters') +parser.add_argument('-r', '--shadow-rotation', type=float, default=0, + help='Arm shadow rotation in degrees,(e.g. 45 degress for some ALMA antennas)') +parser.add_argument('-z', '--first-zscale', type=float, nargs=2, default=[None, None], + help='Z scale for first image (min max)') +parser.add_argument('-y', '--second-zscale', type=float, nargs=2, default=[None, None], + help='Z scale for second image (min max)') +parser.add_argument('-f', '--fits', action='store_true', + default=False, help='Save products as FITS images') + +args = parser.parse_args() + +def get_axis_from_header(header, iaxis): + n_elem = header[f'NAXIS{iaxis}'] + ref = header[f'CRPIX{iaxis}'] + val = header[f'CRVAL{iaxis}'] + inc = header[f'CDELT{iaxis}'] + axis = np.ndarray((n_elem)) + for i_elem in range(n_elem): + axis[i_elem] = val+(ref-i_elem)*inc + return axis + + +def put_axis_in_header(axis, unit, iaxis, header): + n_elem = len(axis) + ref = n_elem//2 + val = axis[ref] + inc = axis[1]-axis[0] + header[f'NAXIS{iaxis}'] = n_elem + header[f'CRPIX{iaxis}'] = ref + header[f'CRVAL{iaxis}'] = val + header[f'CDELT{iaxis}'] = inc + header[f'CUNIT{iaxis}'] = unit + + +def create_fits(header, data, filename): + hdu = fits.PrimaryHDU(data) + for key, value in header.items(): + if isinstance(value, str): + if '/' in value: + wrds = value.split('/') + hdu.header.set(key, wrds[0], wrds[1]) + else: + hdu.header.set(key, value) + hdu.header.set('ORIGIN', f'Image comparison code') + hdu.header.set('DATE', datetime.datetime.now().strftime('%b %d %Y, %H:%M:%S')) + hdu.writeto(filename, overwrite=True) + + +def vers_comp_recursive(ref, cur): + if ref[0] > cur[0]: + return -1 + elif ref[0] < cur[0]: + return 1 + else: + if len(ref) == 1: + return 0 + else: + return vers_comp_recursive(ref[1:], cur[1:]) + + +def test_version(reference, current): + ref_num = get_numbers_from_version(reference) + cur_num = get_numbers_from_version(current) + return vers_comp_recursive(ref_num, cur_num) + + +def get_numbers_from_version(version): + numbers = version[1:].split('.') + revision = int(numbers[0]) + major = int(numbers[1]) + minor = int(numbers[2]) + return [revision, major, minor] + + +class image: + + def __init__(self, filename, noise_level, inlim, oulim, no_division, arm_width, arm_angle): + self.no_division = no_division + self.inlim = inlim + self.oulim = oulim + self.noise_level = noise_level + self.arm_width = arm_width + self.arm_angle = arm_angle*np.pi/180 + + opn_fits = fits.open(filename) + self.data = opn_fits[0].data[0,0,:,:] + self._get_info_from_header(opn_fits[0].header) + opn_fits.close() + + self._mask_image() + + self.division = None + self.residuals = None + self.resampled = False + self.res_mean = None + self.res_rms = None + self.rootname = '.'.join(filename.split('.')[:-1]) + + + def _astrohack_specific_init(self, header): + self.x_unit = header['CUNIT1'] + self.y_unit = header['CUNIT2'] + version = header['ORIGIN'].split()[1][:-1] + if test_version('v0.4.1', version) <= 0: + # This treatment is only necessary before v0.4.2 + wavelength = header['WAVELENG'] + self.x_axis *= wavelength + self.y_axis *= wavelength + self.data = np.fliplr(np.flipud(self.data)) + if test_version('v0.5.0', version) >= 0: + print('estoy aqui?') + wavelength = header['WAVELENG'] + self.x_axis /= wavelength + self.y_axis /= wavelength + self.data = np.fliplr(np.flipud(self.data)) + + + def _get_info_from_header(self, header): + self.header = header + self.x_axis = get_axis_from_header(header, 1) + self.y_axis = get_axis_from_header(header, 2) + self.unit = header['BUNIT'] + + if 'AIPS' in header['ORIGIN']: + self.x_unit = 'm' + self.y_unit = 'm' + elif 'Astrohack' in header['ORIGIN']: + self._astrohack_specific_init(header) + else: + raise Exception(f'Unrecognized origin:\n{header["origin"]}') + + + def _mask_image(self): + x_mesh, y_mesh = np.meshgrid(self.x_axis, self.y_axis) + self.radius = np.sqrt(x_mesh**2 + y_mesh**2) + mask = np.where(self.radius > self.oulim, np.nan, 1.0) + mask = np.where(self.radius < self.inlim, np.nan, mask) + + # Arm masking + if self.arm_angle%np.pi == 0: + mask = np.where(np.abs(x_mesh) < self.arm_width/2., np.nan, mask) + mask = np.where(np.abs(y_mesh) < self.arm_width/2., np.nan, mask) + else: + # first shadow + coeff = np.tan(self.arm_angle%np.pi) + distance = np.abs((coeff*x_mesh-y_mesh)/np.sqrt(coeff**2+1)) + mask = np.where(distance < self.arm_width/2., np.nan, mask) + # second shadow + coeff = np.tan(self.arm_angle%np.pi+np.pi/2) + distance = np.abs((coeff*x_mesh-y_mesh)/np.sqrt(coeff**2+1)) + mask = np.where(distance < self.arm_width/2., np.nan, mask) + + + if self.no_division: + self.noise = None + self.rms = None + else: + self.noise, self.rms = self._noise_filter() + mask = np.where(self.data10*rude_factor, np.nan, division) + self.factor = np.nanmedian(self.division) + + + def _compute_residuals(self, image2, blanking=100): + percent = 100*(self.masked-self.factor*image2.masked)/self.masked + percent = np.where(np.abs(percent)>blanking, np.nan, percent) + self.residuals = percent + self.res_mean = np.nanmean(percent) + self.res_rms = np.sqrt(np.nanmean(percent**2)) + + + def _plot_map(self, data, title, label, filename, cmap, extent, zscale): + fig, ax = plt.subplots(1, 1, figsize=[10,8]) + cmap = get_proper_color_map(cmap) + if zscale == 'symmetrical': + scale = max(np.abs(np.nanmin(data)), np.abs(np.nanmax(data))) + vmin, vmax = -scale, scale + else: + vmin, vmax = zscale + if vmin == 'None' or vmin is None: + vmin = np.nanmin(data) + if vmax == 'None' or vmax is None: + vmax = np.nanmax(data) + + im = ax.imshow(data, cmap=cmap, interpolation="nearest", extent=extent, + vmin=vmin, vmax=vmax,) + well_positioned_colorbar(ax, fig, im, label, location='right', size='5%', pad=0.05) + ax.set_xlabel(f"X axis [{self.x_unit}]") + ax.set_ylabel(f"Y axis [{self.y_unit}]") + close_figure(fig, title, filename, 300, False) + + + def plot(self, plot_noise, cmap, z_scale): + extent = [np.min(self.x_axis), np.max(self.x_axis), np.min(self.y_axis), np.max(self.y_axis)] + png = '.png' + if self.resampled: + title = f'Resampled {self.rootname}' + filename = f'{self.rootname}.resampled' + else: + title = f'{self.rootname}' + filename = f'{self.rootname}' + + if self.no_division: + zlabel = f'?What is type? [{self.unit}]' + else: + zlabel = f'Amplitude [{self.unit}]' + + self._plot_map(self.masked, title, zlabel, filename+png, cmap, extent, z_scale) + self._plot_map(self.mask, f'Mask used for {self.rootname}', 'Mask value', + f'{self.rootname}.mask{png}', cmap, extent, [0, 1]) + + if self.division is not None: + self._plot_map(self.division, f'Division map, mean factor:{self.factor:.3}', + 'Divided value [ ]', + f'{filename}.division{png}', cmap, extent, [None, None]) + + if self.residuals is not None: + self._plot_map(self.residuals, + f'Residual map, mean residual: {self.res_mean:.3}%, residual RMS: ' + f'{self.res_rms:.3}%', + 'Residuals [%]', + f'{filename}.residual{png}', 'RdBu_r', extent, 'symmetrical') + + if plot_noise: + if self.noise is not None: + self._plot_map(self.noise, + f'Noise like component for {self.rootname}', zlabel, + f'{filename}.noise{png}', cmap, extent, [None, None]) + + + def to_fits(self): + fits = '.fits' + header = self.header.copy() + put_axis_in_header(self.x_axis, self.x_unit, 1, header) + put_axis_in_header(self.y_axis, self.y_unit, 2, header) + + if self.resampled: + filename = f'comp_{self.rootname}.resampled' + else: + filename = f'comp_{self.rootname}' + + create_fits(header, self.masked, f'{filename}.masked{fits}') + create_fits(header, self.mask, f'{filename}.mask{fits}') + + if self.division is not None: + create_fits(header, self.division, f'{filename}.division{fits}') + + if self.residuals is not None: + create_fits(header, self.residuals, f'{filename}.residual{fits}') + + + if self.noise is not None: + create_fits(header, self.noise, f'{filename}.noise{fits}') + + + def print_stats(self): + print(80*'*') + print() + print(f'Mean scaling factor = {self.factor:.3}') + print(f'Mean Residual = {self.res_mean:.3}%') + print(f'Residuals RMS = {self.res_rms:.3}%') + print() + + +# instatiation +first_image = image(args.first, args.noise_clip, args.blocage, + args.diameter/2, args.no_division, args.shadow_width, + args.shadow_rotation) +second_image = image(args.second, args.noise_clip, args.blocage, + args.diameter/2, args.no_division, args.shadow_width, + args.shadow_rotation) + +# Data manipulation +second_image.resample(first_image) +first_image.make_comparison(second_image) + +# Plotting +first_image.plot(args.noise_map, args.colormap, args.first_zscale) +second_image.plot(args.noise_map, args.colormap, args.second_zscale) + +if args.fits: + first_image.to_fits() + second_image.to_fits() + +if not args.quiet: + first_image.print_stats() + + + + + From 4a5a139d3e6f1d9a8945639e4dc99d7e2d24dfc2 Mon Sep 17 00:00:00 2001 From: Victor de Souza Magalhaes Date: Tue, 18 Feb 2025 14:53:32 -0700 Subject: [PATCH 07/53] Added missing docstrings in fits.py --- src/astrohack/utils/fits.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/astrohack/utils/fits.py b/src/astrohack/utils/fits.py index d271f4a8..fc3b726d 100644 --- a/src/astrohack/utils/fits.py +++ b/src/astrohack/utils/fits.py @@ -103,6 +103,8 @@ def axis_to_fits_header(header: dict, axis, iaxis, axistype, unit, iswcs=True): axis: The axis to be described in the header iaxis: The position of the axis in the data axistype: Axis type to be displayed in the fits header + unit: Axis unit + iswcs: Is the axis a part of World Coordinate System for the image? Returns: The augmented header From 240e2d893dd80d5d6aca9cf30e54915f9ff59a96 Mon Sep 17 00:00:00 2001 From: Victor de Souza Magalhaes Date: Tue, 18 Feb 2025 14:58:48 -0700 Subject: [PATCH 08/53] Removed obviously redundant Ad hoc solutions from image_compare_tool.py --- src/astrohack/core/image_compare_tool.py | 151 +++++------------------ 1 file changed, 33 insertions(+), 118 deletions(-) diff --git a/src/astrohack/core/image_compare_tool.py b/src/astrohack/core/image_compare_tool.py index 17ff0c0e..803eddee 100644 --- a/src/astrohack/core/image_compare_tool.py +++ b/src/astrohack/core/image_compare_tool.py @@ -2,107 +2,54 @@ import numpy as np from scipy.interpolate import griddata from matplotlib import pyplot as plt -from mpl_toolkits.axes_grid1 import make_axes_locatable -import argparse from astrohack.visualization.plot_tools import well_positioned_colorbar from astrohack.visualization.plot_tools import close_figure, get_proper_color_map +from astrohack.utils.fits import read_fits, axis_to_fits_header, write_fits import datetime -parser = argparse.ArgumentParser(description="Compare aperture FTIS maps produced by AIPS and AstroHACK\n" - "Beam maps are not supported") -parser.add_argument('first', type=str, help='first aperture FITS (usually AIPS)') -parser.add_argument('second', type=str, - help='Second aperture FITS (usually AstroHACK)') -parser.add_argument('-n', '--noise-map', action='store_true', - default=False, help='Save noise images') -parser.add_argument('-c', '--noise-clip', type=float, default=1, - help='Noise clipping level (in sigmas)') -parser.add_argument('-q', '--quiet', action='store_true', default=False, - help='Do not print value to screen') -parser.add_argument('-d', '--diameter', type=float, default=25.0, - help='Dish diameter') -parser.add_argument('-b', '--blocage', type=float, default=2.0, - help='Dish inner blocage radius') -parser.add_argument('-m', '--colormap', type=str, default='viridis', - help='Colormap for non residual maps') -parser.add_argument('-p', '--no-division', action='store_true', default=False, - help='Do not perform division, i.e. factor assumed to be 1') -parser.add_argument('-s', '--shadow-width', type=float, default=1.5, - help='Arm shadow width in meters') -parser.add_argument('-r', '--shadow-rotation', type=float, default=0, - help='Arm shadow rotation in degrees,(e.g. 45 degress for some ALMA antennas)') -parser.add_argument('-z', '--first-zscale', type=float, nargs=2, default=[None, None], - help='Z scale for first image (min max)') -parser.add_argument('-y', '--second-zscale', type=float, nargs=2, default=[None, None], - help='Z scale for second image (min max)') -parser.add_argument('-f', '--fits', action='store_true', - default=False, help='Save products as FITS images') - -args = parser.parse_args() +# parser = argparse.ArgumentParser(description="Compare aperture FTIS maps produced by AIPS and AstroHACK\n" +# "Beam maps are not supported") +# parser.add_argument('first', type=str, help='first aperture FITS (usually AIPS)') +# parser.add_argument('second', type=str, +# help='Second aperture FITS (usually AstroHACK)') +# parser.add_argument('-n', '--noise-map', action='store_true', +# default=False, help='Save noise images') +# parser.add_argument('-c', '--noise-clip', type=float, default=1, +# help='Noise clipping level (in sigmas)') +# parser.add_argument('-q', '--quiet', action='store_true', default=False, +# help='Do not print value to screen') +# parser.add_argument('-d', '--diameter', type=float, default=25.0, +# help='Dish diameter') +# parser.add_argument('-b', '--blocage', type=float, default=2.0, +# help='Dish inner blocage radius') +# parser.add_argument('-m', '--colormap', type=str, default='viridis', +# help='Colormap for non residual maps') +# parser.add_argument('-p', '--no-division', action='store_true', default=False, +# help='Do not perform division, i.e. factor assumed to be 1') +# parser.add_argument('-s', '--shadow-width', type=float, default=1.5, +# help='Arm shadow width in meters') +# parser.add_argument('-r', '--shadow-rotation', type=float, default=0, +# help='Arm shadow rotation in degrees,(e.g. 45 degress for some ALMA antennas)') +# parser.add_argument('-z', '--first-zscale', type=float, nargs=2, default=[None, None], +# help='Z scale for first image (min max)') +# parser.add_argument('-y', '--second-zscale', type=float, nargs=2, default=[None, None], +# help='Z scale for second image (min max)') +# parser.add_argument('-f', '--fits', action='store_true', +# default=False, help='Save products as FITS images') +# +# args = parser.parse_args() def get_axis_from_header(header, iaxis): n_elem = header[f'NAXIS{iaxis}'] ref = header[f'CRPIX{iaxis}'] val = header[f'CRVAL{iaxis}'] inc = header[f'CDELT{iaxis}'] - axis = np.ndarray((n_elem)) + axis = np.ndarray(n_elem) for i_elem in range(n_elem): axis[i_elem] = val+(ref-i_elem)*inc return axis -def put_axis_in_header(axis, unit, iaxis, header): - n_elem = len(axis) - ref = n_elem//2 - val = axis[ref] - inc = axis[1]-axis[0] - header[f'NAXIS{iaxis}'] = n_elem - header[f'CRPIX{iaxis}'] = ref - header[f'CRVAL{iaxis}'] = val - header[f'CDELT{iaxis}'] = inc - header[f'CUNIT{iaxis}'] = unit - - -def create_fits(header, data, filename): - hdu = fits.PrimaryHDU(data) - for key, value in header.items(): - if isinstance(value, str): - if '/' in value: - wrds = value.split('/') - hdu.header.set(key, wrds[0], wrds[1]) - else: - hdu.header.set(key, value) - hdu.header.set('ORIGIN', f'Image comparison code') - hdu.header.set('DATE', datetime.datetime.now().strftime('%b %d %Y, %H:%M:%S')) - hdu.writeto(filename, overwrite=True) - - -def vers_comp_recursive(ref, cur): - if ref[0] > cur[0]: - return -1 - elif ref[0] < cur[0]: - return 1 - else: - if len(ref) == 1: - return 0 - else: - return vers_comp_recursive(ref[1:], cur[1:]) - - -def test_version(reference, current): - ref_num = get_numbers_from_version(reference) - cur_num = get_numbers_from_version(current) - return vers_comp_recursive(ref_num, cur_num) - - -def get_numbers_from_version(version): - numbers = version[1:].split('.') - revision = int(numbers[0]) - major = int(numbers[1]) - minor = int(numbers[2]) - return [revision, major, minor] - - class image: def __init__(self, filename, noise_level, inlim, oulim, no_division, arm_width, arm_angle): @@ -161,38 +108,6 @@ def _get_info_from_header(self, header): raise Exception(f'Unrecognized origin:\n{header["origin"]}') - def _mask_image(self): - x_mesh, y_mesh = np.meshgrid(self.x_axis, self.y_axis) - self.radius = np.sqrt(x_mesh**2 + y_mesh**2) - mask = np.where(self.radius > self.oulim, np.nan, 1.0) - mask = np.where(self.radius < self.inlim, np.nan, mask) - - # Arm masking - if self.arm_angle%np.pi == 0: - mask = np.where(np.abs(x_mesh) < self.arm_width/2., np.nan, mask) - mask = np.where(np.abs(y_mesh) < self.arm_width/2., np.nan, mask) - else: - # first shadow - coeff = np.tan(self.arm_angle%np.pi) - distance = np.abs((coeff*x_mesh-y_mesh)/np.sqrt(coeff**2+1)) - mask = np.where(distance < self.arm_width/2., np.nan, mask) - # second shadow - coeff = np.tan(self.arm_angle%np.pi+np.pi/2) - distance = np.abs((coeff*x_mesh-y_mesh)/np.sqrt(coeff**2+1)) - mask = np.where(distance < self.arm_width/2., np.nan, mask) - - - if self.no_division: - self.noise = None - self.rms = None - else: - self.noise, self.rms = self._noise_filter() - mask = np.where(self.data Date: Thu, 20 Feb 2025 14:43:59 -0700 Subject: [PATCH 09/53] Renamed some FITS routines for clearer semantics. --- src/astrohack/antenna/antenna_surface.py | 8 ++++---- src/astrohack/core/image_compare_tool.py | 2 +- src/astrohack/utils/fits.py | 6 +++--- src/astrohack/visualization/fits.py | 22 +++++++++++----------- 4 files changed, 19 insertions(+), 19 deletions(-) diff --git a/src/astrohack/antenna/antenna_surface.py b/src/astrohack/antenna/antenna_surface.py index f1aac746..acd2b7fb 100644 --- a/src/astrohack/antenna/antenna_surface.py +++ b/src/astrohack/antenna/antenna_surface.py @@ -14,7 +14,7 @@ from astrohack.visualization.plot_tools import well_positioned_colorbar, create_figure_and_axes, close_figure, \ get_proper_color_map -from astrohack.utils.fits import write_fits, resolution_to_fits_header, axis_to_fits_header +from astrohack.utils.fits import write_fits, put_resolution_in_fits_header, put_axis_in_fits_header lnbr = "\n" SUPPORTED_POL_STATES = ['I', 'RR', 'LL', 'XX', 'YY'] @@ -865,9 +865,9 @@ def export_to_fits(self, basename): 'WAVELENG': self.wavelength, 'FREQUENC': clight / self.wavelength, } - head = axis_to_fits_header(head, self.u_axis, 1, 'X----LIN', 'm') - head = axis_to_fits_header(head, self.v_axis, 2, 'Y----LIN', 'm') - head = resolution_to_fits_header(head, self.resolution) + head = put_axis_in_fits_header(head, self.u_axis, 1, 'X----LIN', 'm') + head = put_axis_in_fits_header(head, self.v_axis, 2, 'Y----LIN', 'm') + head = put_resolution_in_fits_header(head, self.resolution) write_fits(head, 'Amplitude', self.amplitude, add_prefix(basename, 'amplitude') + '.fits', self.amp_unit, 'panel') diff --git a/src/astrohack/core/image_compare_tool.py b/src/astrohack/core/image_compare_tool.py index 803eddee..88f98eb3 100644 --- a/src/astrohack/core/image_compare_tool.py +++ b/src/astrohack/core/image_compare_tool.py @@ -4,7 +4,7 @@ from matplotlib import pyplot as plt from astrohack.visualization.plot_tools import well_positioned_colorbar from astrohack.visualization.plot_tools import close_figure, get_proper_color_map -from astrohack.utils.fits import read_fits, axis_to_fits_header, write_fits +from astrohack.utils.fits import read_fits, put_axis_in_fits_header, write_fits import datetime # parser = argparse.ArgumentParser(description="Compare aperture FTIS maps produced by AIPS and AstroHACK\n" diff --git a/src/astrohack/utils/fits.py b/src/astrohack/utils/fits.py index fc3b726d..05299929 100644 --- a/src/astrohack/utils/fits.py +++ b/src/astrohack/utils/fits.py @@ -73,7 +73,7 @@ def _reorder_axes_for_fits(data: np.ndarray): return np.flipud(data) -def resolution_to_fits_header(header, resolution): +def put_resolution_in_fits_header(header, resolution): """ Adds resolution information to standard header keywords: BMAJ, BMIN and BPA Args: @@ -95,7 +95,7 @@ def resolution_to_fits_header(header, resolution): return header -def axis_to_fits_header(header: dict, axis, iaxis, axistype, unit, iswcs=True): +def put_axis_in_fits_header(header: dict, axis, iaxis, axistype, unit, iswcs=True): """ Process an axis to create a FITS compatible linear axis description Args: @@ -144,7 +144,7 @@ def axis_to_fits_header(header: dict, axis, iaxis, axistype, unit, iswcs=True): return outheader -def stokes_axis_to_fits_header(header, iaxis): +def put_stokes_axis_in_fits_header(header, iaxis): """ Inserts a dedicated stokes axis in the header at iaxis Args: diff --git a/src/astrohack/visualization/fits.py b/src/astrohack/visualization/fits.py index 1db92166..c1dcd4b8 100644 --- a/src/astrohack/visualization/fits.py +++ b/src/astrohack/visualization/fits.py @@ -2,7 +2,7 @@ from toolviper.utils import logger as logger from astrohack.antenna import Telescope, AntennaSurface from astrohack.utils import clight, convert_unit, add_prefix -from astrohack.utils.fits import axis_to_fits_header, stokes_axis_to_fits_header, write_fits, resolution_to_fits_header +from astrohack.utils.fits import put_axis_in_fits_header, put_stokes_axis_in_fits_header, write_fits, put_resolution_in_fits_header def export_to_fits_panel_chunk(parm_dict): @@ -85,11 +85,11 @@ def export_to_fits_holog_chunk(parm_dict): if ntime != 1: raise Exception("Data with multiple times not supported for FITS export") - base_header = axis_to_fits_header(base_header, input_xds.chan.values, 3, 'Frequency', 'Hz') - base_header = stokes_axis_to_fits_header(base_header, 4) + base_header = put_axis_in_fits_header(base_header, input_xds.chan.values, 3, 'Frequency', 'Hz') + base_header = put_stokes_axis_in_fits_header(base_header, 4) rad_to_deg = convert_unit('rad', 'deg', 'trigonometric') - beam_header = axis_to_fits_header(base_header, -input_xds.l.values * rad_to_deg, 1, 'RA---SIN', 'deg') - beam_header = axis_to_fits_header(beam_header, input_xds.m.values * rad_to_deg, 2, 'DEC--SIN', 'deg') + beam_header = put_axis_in_fits_header(base_header, -input_xds.l.values * rad_to_deg, 1, 'RA---SIN', 'deg') + beam_header = put_axis_in_fits_header(beam_header, input_xds.m.values * rad_to_deg, 2, 'DEC--SIN', 'deg') beam_header['RADESYSA'] = 'FK5' beam = input_xds['BEAM'].values if parm_dict['complex_split'] == 'cartesian': @@ -103,9 +103,9 @@ def export_to_fits_holog_chunk(parm_dict): write_fits(beam_header, 'Complex beam phase', np.angle(beam), add_prefix(basename, 'beam_phase') + '.fits', 'Radians', 'image') wavelength = clight / input_xds.chan.values[0] - aperture_header = axis_to_fits_header(base_header, input_xds.u.values * wavelength, 1, 'X----LIN', 'm') - aperture_header = axis_to_fits_header(aperture_header, input_xds.u.values * wavelength, 2, 'Y----LIN', 'm') - aperture_header = resolution_to_fits_header(aperture_header, aperture_resolution) + aperture_header = put_axis_in_fits_header(base_header, input_xds.u.values * wavelength, 1, 'X----LIN', 'm') + aperture_header = put_axis_in_fits_header(aperture_header, input_xds.u.values * wavelength, 2, 'Y----LIN', 'm') + aperture_header = put_resolution_in_fits_header(aperture_header, aperture_resolution) aperture = input_xds['APERTURE'].values if parm_dict['complex_split'] == 'cartesian': write_fits(aperture_header, 'Complex aperture real part', aperture.real, @@ -118,9 +118,9 @@ def export_to_fits_holog_chunk(parm_dict): write_fits(aperture_header, 'Complex aperture phase', np.angle(aperture), add_prefix(basename, 'aperture_phase') + '.fits', 'rad', 'image') - phase_amp_header = axis_to_fits_header(base_header, input_xds.u_prime.values * wavelength, 1, 'X----LIN', 'm') - phase_amp_header = axis_to_fits_header(phase_amp_header, input_xds.v_prime.values * wavelength, 2, 'Y----LIN', 'm') - phase_amp_header = resolution_to_fits_header(phase_amp_header, aperture_resolution) + phase_amp_header = put_axis_in_fits_header(base_header, input_xds.u_prime.values * wavelength, 1, 'X----LIN', 'm') + phase_amp_header = put_axis_in_fits_header(phase_amp_header, input_xds.v_prime.values * wavelength, 2, 'Y----LIN', 'm') + phase_amp_header = put_resolution_in_fits_header(phase_amp_header, aperture_resolution) write_fits(phase_amp_header, 'Cropped aperture corrected phase', input_xds['CORRECTED_PHASE'].values, add_prefix(basename, 'corrected_phase') + '.fits', 'rad', 'image') return From a617bd41b959e9ef60cc270097077e23d65bb756 Mon Sep 17 00:00:00 2001 From: Victor de Souza magalhaes Date: Thu, 20 Feb 2025 14:53:00 -0700 Subject: [PATCH 10/53] Moved get_axis_from_fits_header to utils/fits.py --- src/astrohack/core/image_compare_tool.py | 10 ---------- src/astrohack/utils/fits.py | 21 +++++++++++++++++++++ 2 files changed, 21 insertions(+), 10 deletions(-) diff --git a/src/astrohack/core/image_compare_tool.py b/src/astrohack/core/image_compare_tool.py index 88f98eb3..4cff49f2 100644 --- a/src/astrohack/core/image_compare_tool.py +++ b/src/astrohack/core/image_compare_tool.py @@ -39,16 +39,6 @@ # # args = parser.parse_args() -def get_axis_from_header(header, iaxis): - n_elem = header[f'NAXIS{iaxis}'] - ref = header[f'CRPIX{iaxis}'] - val = header[f'CRVAL{iaxis}'] - inc = header[f'CDELT{iaxis}'] - axis = np.ndarray(n_elem) - for i_elem in range(n_elem): - axis[i_elem] = val+(ref-i_elem)*inc - return axis - class image: diff --git a/src/astrohack/utils/fits.py b/src/astrohack/utils/fits.py index 05299929..17682fc3 100644 --- a/src/astrohack/utils/fits.py +++ b/src/astrohack/utils/fits.py @@ -35,6 +35,27 @@ def read_fits(filename): return head, data +def get_axis_from_fits_header(header, iaxis): + """ + Pull axis information from FITS file and store it in a numpy array, ignores rotation in axes. + Args: + header: FITS header + iaxis: Which axis is to be fetched from the header. + + Returns: + numpy array representation of axis, axis type and axis unit + """ + n_elem = header[f'NAXIS{iaxis}'] + ref = header[f'CRPIX{iaxis}'] + inc = header[f'CDELT{iaxis}'] + val = header[f'CRVAL{iaxis}'] + inc # This makes this routine symmetrical to the put routine. + axis = np.arange(n_elem) + axis = val + (ref-axis)*inc + axis_unit = header[f'CUNIT{iaxis}'] + axis_type = header[f'CTYPE{iaxis}'] + return axis, axis_type, axis_unit + + def write_fits(header, imagetype, data, filename, unit, origin): """ Write a dictionary and a dataset to a FITS file From c0921615a94295434864033b69bf4e944bf660a4 Mon Sep 17 00:00:00 2001 From: Victor de Souza magalhaes Date: Thu, 20 Feb 2025 16:05:39 -0700 Subject: [PATCH 11/53] Added a simple tool to compare axes --- src/astrohack/utils/algorithms.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/astrohack/utils/algorithms.py b/src/astrohack/utils/algorithms.py index bf9e2f47..957a36d1 100644 --- a/src/astrohack/utils/algorithms.py +++ b/src/astrohack/utils/algorithms.py @@ -591,3 +591,9 @@ def _arm_shadow_masking(inmask, x_mesh, y_mesh, radius_mesh, minradius, maxradiu distance = np.abs((coeff*x_mesh-y_mesh)/np.sqrt(coeff**2+1)) oumask = np.where(np.bitwise_and(distance < width/2., radial_mask), False, oumask) return oumask + + +def are_axes_equal(axis_a, axis_b): + if axis_a.shape[0] != axis_b.shape[0]: + return False + return np.all(axis_a == axis_b) From 8c0c8ead074cab55c9403139eade2deb8107cf89 Mon Sep 17 00:00:00 2001 From: Victor de Souza magalhaes Date: Thu, 20 Feb 2025 16:06:16 -0700 Subject: [PATCH 12/53] Implemented comparison methods and plotting --- src/astrohack/core/image_compare_tool.py | 326 +++++++++-------------- 1 file changed, 132 insertions(+), 194 deletions(-) diff --git a/src/astrohack/core/image_compare_tool.py b/src/astrohack/core/image_compare_tool.py index 4cff49f2..0d807bb7 100644 --- a/src/astrohack/core/image_compare_tool.py +++ b/src/astrohack/core/image_compare_tool.py @@ -1,153 +1,127 @@ -from astropy.io import fits import numpy as np from scipy.interpolate import griddata from matplotlib import pyplot as plt + +from astrohack import compute_extent, data_statistics, statistics_to_text +from astrohack.utils import are_axes_equal +from astrohack.utils.algorithms import create_aperture_mask from astrohack.visualization.plot_tools import well_positioned_colorbar from astrohack.visualization.plot_tools import close_figure, get_proper_color_map -from astrohack.utils.fits import read_fits, put_axis_in_fits_header, write_fits +from astrohack.utils.fits import read_fits, put_axis_in_fits_header, write_fits, get_axis_from_fits_header import datetime -# parser = argparse.ArgumentParser(description="Compare aperture FTIS maps produced by AIPS and AstroHACK\n" -# "Beam maps are not supported") -# parser.add_argument('first', type=str, help='first aperture FITS (usually AIPS)') -# parser.add_argument('second', type=str, -# help='Second aperture FITS (usually AstroHACK)') -# parser.add_argument('-n', '--noise-map', action='store_true', -# default=False, help='Save noise images') -# parser.add_argument('-c', '--noise-clip', type=float, default=1, -# help='Noise clipping level (in sigmas)') -# parser.add_argument('-q', '--quiet', action='store_true', default=False, -# help='Do not print value to screen') -# parser.add_argument('-d', '--diameter', type=float, default=25.0, -# help='Dish diameter') -# parser.add_argument('-b', '--blocage', type=float, default=2.0, -# help='Dish inner blocage radius') -# parser.add_argument('-m', '--colormap', type=str, default='viridis', -# help='Colormap for non residual maps') -# parser.add_argument('-p', '--no-division', action='store_true', default=False, -# help='Do not perform division, i.e. factor assumed to be 1') -# parser.add_argument('-s', '--shadow-width', type=float, default=1.5, -# help='Arm shadow width in meters') -# parser.add_argument('-r', '--shadow-rotation', type=float, default=0, -# help='Arm shadow rotation in degrees,(e.g. 45 degress for some ALMA antennas)') -# parser.add_argument('-z', '--first-zscale', type=float, nargs=2, default=[None, None], -# help='Z scale for first image (min max)') -# parser.add_argument('-y', '--second-zscale', type=float, nargs=2, default=[None, None], -# help='Z scale for second image (min max)') -# parser.add_argument('-f', '--fits', action='store_true', -# default=False, help='Save products as FITS images') -# -# args = parser.parse_args() - - -class image: - - def __init__(self, filename, noise_level, inlim, oulim, no_division, arm_width, arm_angle): - self.no_division = no_division - self.inlim = inlim - self.oulim = oulim - self.noise_level = noise_level - self.arm_width = arm_width - self.arm_angle = arm_angle*np.pi/180 - - opn_fits = fits.open(filename) - self.data = opn_fits[0].data[0,0,:,:] - self._get_info_from_header(opn_fits[0].header) - opn_fits.close() - - self._mask_image() - - self.division = None - self.residuals = None - self.resampled = False - self.res_mean = None - self.res_rms = None - self.rootname = '.'.join(filename.split('.')[:-1]) +def test_image(fits_image): + if isinstance(fits_image, FITSImage): + pass + else: + raise TypeError('Reference image is not a FITSImage object') + + +class FITSImage: - def _astrohack_specific_init(self, header): - self.x_unit = header['CUNIT1'] - self.y_unit = header['CUNIT2'] - version = header['ORIGIN'].split()[1][:-1] - if test_version('v0.4.1', version) <= 0: - # This treatment is only necessary before v0.4.2 - wavelength = header['WAVELENG'] - self.x_axis *= wavelength - self.y_axis *= wavelength - self.data = np.fliplr(np.flipud(self.data)) - if test_version('v0.5.0', version) >= 0: - print('estoy aqui?') - wavelength = header['WAVELENG'] - self.x_axis /= wavelength - self.y_axis /= wavelength - self.data = np.fliplr(np.flipud(self.data)) - - - def _get_info_from_header(self, header): - self.header = header - self.x_axis = get_axis_from_header(header, 1) - self.y_axis = get_axis_from_header(header, 2) - self.unit = header['BUNIT'] - - if 'AIPS' in header['ORIGIN']: + def __init__(self, filename, telescope_obj): + self.telescope = telescope_obj + self.header, data = read_fits(filename) + self.data = data[0, 0, :, :] + + self.x_axis, _, self.x_unit = get_axis_from_fits_header(self.header, 1) + self.y_axis, _, self.y_unit = get_axis_from_fits_header(self.header, 2) + + self.unit = self.header['BUNIT'] + + if 'AIPS' in self.header['ORIGIN']: self.x_unit = 'm' self.y_unit = 'm' - elif 'Astrohack' in header['ORIGIN']: - self._astrohack_specific_init(header) + elif 'Astrohack' in self.header['ORIGIN']: + self.data = self.data.T else: - raise Exception(f'Unrecognized origin:\n{header["origin"]}') + raise Exception(f'Unrecognized origin:\n{self.header["origin"]}') + self.base_mask = create_aperture_mask(self.x_axis, self.y_axis, self.telescope.inlim, self.telescope.oulim, + arm_width=self.telescope.arm_shadow_width, + arm_angle=self.telescope.arm_shadow_rotation) + self.rootname = '.'.join(filename.split('.')[:-1]) + self.factor = 1.0 + self.residuals = None + self.residuals_percent = None + self.divided_image = None + self.reference_name = None def resample(self, ref_image): - x_mesh_orig, y_mesh_orig = np.meshgrid(self.x_axis, self.y_axis) - x_mesh_dest, y_mesh_dest = np.meshgrid(ref_image.x_axis, ref_image.y_axis) - resamp = griddata((x_mesh_orig.ravel(), y_mesh_orig.ravel()), - self.data.ravel(), + test_image(ref_image) + x_mesh_orig, y_mesh_orig = np.meshgrid(self.x_axis, self.y_axis, indexing='ij') + x_mesh_dest, y_mesh_dest = np.meshgrid(ref_image.x_axis, ref_image.y_axis, indexing='ij') + resamp = griddata((x_mesh_orig.ravel(), y_mesh_orig.ravel()), self.data.ravel(), (x_mesh_dest.ravel(), y_mesh_dest.ravel()), method='linear') size = ref_image.x_axis.shape[0], ref_image.y_axis.shape[0] self.x_axis = ref_image.x_axis self.y_axis = ref_image.y_axis self.data = resamp.reshape(size) - self._mask_image() - self.resampled = True - - - def _noise_filter(self): - noise = np.where(self.radius < self.oulim, np.nan, self.data) - noise = np.where(self.radius < self.inlim, self.data, noise) - noiserms = np.sqrt(np.nanmean(noise**2)) - return noise, noiserms - - - def make_comparison(self, image2): - if self.no_division: - self.factor = 1.0 - self.division = None - else: - self._compute_divide(image2) - - self._compute_residuals(image2) - - - def _compute_divide(self, image2): - division = self.masked/image2.masked - rude_factor = np.abs(np.nanmean(division)) - self.division = np.where(np.abs(division)>10*rude_factor, np.nan, division) - self.factor = np.nanmedian(self.division) - - - def _compute_residuals(self, image2, blanking=100): - percent = 100*(self.masked-self.factor*image2.masked)/self.masked - percent = np.where(np.abs(percent)>blanking, np.nan, percent) - self.residuals = percent - self.res_mean = np.nanmean(percent) - self.res_rms = np.sqrt(np.nanmean(percent**2)) - - - def _plot_map(self, data, title, label, filename, cmap, extent, zscale): - fig, ax = plt.subplots(1, 1, figsize=[10,8]) - cmap = get_proper_color_map(cmap) + self.base_mask = create_aperture_mask(self.x_axis, self.y_axis, self.telescope.inlim, self.telescope.oulim, + arm_width=self.telescope.arm_shadow_width, + arm_angle=self.telescope.arm_shadow_rotation) + + def compare_difference(self, ref_image): + test_image(ref_image) + if not self.image_has_same_sampling(ref_image): + self.resample(ref_image) + + self.residuals = ref_image.data - (self.data * self.factor) + self.residuals_percent = 100 * self.residuals/ref_image.data + self.reference_name = ref_image.rootname + + def compare_scaled_difference(self, ref_image, rejection=10): + test_image(ref_image) + if not self.image_has_same_sampling(ref_image): + self.resample(ref_image) + simple_division = ref_image.data / self.data + rough_factor = np.nanmean(simple_division[self.base_mask]) + self.divided_image = np.where(np.abs(simple_division) > rejection*rough_factor, np.nan, simple_division) + self.factor = np.nanmedian(self.divided_image) + self.compare_difference(ref_image) + + def image_has_same_sampling(self, ref_image): + test_image(ref_image) + return are_axes_equal(self.x_axis, ref_image.x_axis) and are_axes_equal(self.y_axis, ref_image.y_axis) + + def _mask_array(self, image_array): + return np.where(self.base_mask, image_array, np.nan) + + def plot_results(self, destination, plot_data=False, plot_percentuals=False, plot_divided_image=False, + colormap='viridis', dpi=300, display=False): + + extent = compute_extent(self.x_axis, self.y_axis, 0.0) + cmap = get_proper_color_map(colormap) + base_name = f'{destination}/{self.rootname}' + + if self.residuals is None: + raise Exception("Cannot plot results as they don't exist yet.") + + self._plot_map(self._mask_array(self.residuals), 'Residuals', f'Residuals [{self.unit}]', + f'{base_name}residuals.png', cmap, extent, 'Symmetrical', dpi, display, + add_statistics=True) + + if plot_data: + self._plot_map(self._mask_array(self.residuals), 'Original Data', f'Data [{self.unit}]', + f'{base_name}data.png', cmap, extent, [None, None], dpi, display, + add_statistics=False) + + if plot_percentuals: + self._plot_map(self._mask_array(self.residuals_percent), 'Residuals in %', f'Residuals [%]', + f'{base_name}residuals_percent.png', cmap, extent, 'Symmetrical', dpi, display, + add_statistics=True) + + if plot_divided_image: + if self.divided_image is None: + raise Exception("Cannot plot a divided image that does not exist.") + self._plot_map(self._mask_array(self.divided_image), 'Divided image', f'Division [ ]', + f'{base_name}divided.png', cmap, extent, [None, None], dpi, display, + add_statistics=True) + + def _plot_map(self, data, title, zlabel, filename, cmap, extent, zscale, dpi, display, add_statistics=False): + fig, ax = plt.subplots(1, 1, figsize=[10, 8]) if zscale == 'symmetrical': scale = max(np.abs(np.nanmin(data)), np.abs(np.nanmax(data))) vmin, vmax = -scale, scale @@ -160,73 +134,37 @@ def _plot_map(self, data, title, label, filename, cmap, extent, zscale): im = ax.imshow(data, cmap=cmap, interpolation="nearest", extent=extent, vmin=vmin, vmax=vmax,) - well_positioned_colorbar(ax, fig, im, label, location='right', size='5%', pad=0.05) + well_positioned_colorbar(ax, fig, im, zlabel, location='right', size='5%', pad=0.05) ax.set_xlabel(f"X axis [{self.x_unit}]") ax.set_ylabel(f"Y axis [{self.y_unit}]") - close_figure(fig, title, filename, 300, False) - - - def plot(self, plot_noise, cmap, z_scale): - extent = [np.min(self.x_axis), np.max(self.x_axis), np.min(self.y_axis), np.max(self.y_axis)] - png = '.png' - if self.resampled: - title = f'Resampled {self.rootname}' - filename = f'{self.rootname}.resampled' - else: - title = f'{self.rootname}' - filename = f'{self.rootname}' - - if self.no_division: - zlabel = f'?What is type? [{self.unit}]' - else: - zlabel = f'Amplitude [{self.unit}]' - - self._plot_map(self.masked, title, zlabel, filename+png, cmap, extent, z_scale) - self._plot_map(self.mask, f'Mask used for {self.rootname}', 'Mask value', - f'{self.rootname}.mask{png}', cmap, extent, [0, 1]) - - if self.division is not None: - self._plot_map(self.division, f'Division map, mean factor:{self.factor:.3}', - 'Divided value [ ]', - f'{filename}.division{png}', cmap, extent, [None, None]) - - if self.residuals is not None: - self._plot_map(self.residuals, - f'Residual map, mean residual: {self.res_mean:.3}%, residual RMS: ' - f'{self.res_rms:.3}%', - 'Residuals [%]', - f'{filename}.residual{png}', 'RdBu_r', extent, 'symmetrical') - - if plot_noise: - if self.noise is not None: - self._plot_map(self.noise, - f'Noise like component for {self.rootname}', zlabel, - f'{filename}.noise{png}', cmap, extent, [None, None]) - - - def to_fits(self): - fits = '.fits' - header = self.header.copy() - put_axis_in_header(self.x_axis, self.x_unit, 1, header) - put_axis_in_header(self.y_axis, self.y_unit, 2, header) - - if self.resampled: - filename = f'comp_{self.rootname}.resampled' - else: - filename = f'comp_{self.rootname}' - - create_fits(header, self.masked, f'{filename}.masked{fits}') - create_fits(header, self.mask, f'{filename}.mask{fits}') - - if self.division is not None: - create_fits(header, self.division, f'{filename}.division{fits}') - - if self.residuals is not None: - create_fits(header, self.residuals, f'{filename}.residual{fits}') - - - if self.noise is not None: - create_fits(header, self.noise, f'{filename}.noise{fits}') + if add_statistics: + data_stats = data_statistics(data) + ax.set_title(statistics_to_text(data_stats)) + close_figure(fig, title, filename, dpi, display) + + # def to_fits(self): + # fits = '.fits' + # header = self.header.copy() + # put_axis_in_header(self.x_axis, self.x_unit, 1, header) + # put_axis_in_header(self.y_axis, self.y_unit, 2, header) + # + # if self.resampled: + # filename = f'comp_{self.rootname}.resampled' + # else: + # filename = f'comp_{self.rootname}' + # + # create_fits(header, self.masked, f'{filename}.masked{fits}') + # create_fits(header, self.mask, f'{filename}.mask{fits}') + # + # if self.division is not None: + # create_fits(header, self.division, f'{filename}.division{fits}') + # + # if self.residuals is not None: + # create_fits(header, self.residuals, f'{filename}.residual{fits}') + # + # + # if self.noise is not None: + # create_fits(header, self.noise, f'{filename}.noise{fits}') def print_stats(self): From 1a41a04f9af3ee594e7bf5d3aa89a56c976223ff Mon Sep 17 00:00:00 2001 From: Victor de Souza magalhaes Date: Thu, 20 Feb 2025 16:41:58 -0700 Subject: [PATCH 13/53] Evolved FITS reading and writing. --- src/astrohack/utils/fits.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/astrohack/utils/fits.py b/src/astrohack/utils/fits.py index 17682fc3..e6326304 100644 --- a/src/astrohack/utils/fits.py +++ b/src/astrohack/utils/fits.py @@ -21,7 +21,7 @@ def read_fits(filename): """ hdul = fits.open(filename) head = hdul[0].header - data = hdul[0].data[0, 0, :, :] + data = hdul[0].data hdul.close() if head["NAXIS"] != 1: if head["NAXIS"] < 1: @@ -56,7 +56,7 @@ def get_axis_from_fits_header(header, iaxis): return axis, axis_type, axis_unit -def write_fits(header, imagetype, data, filename, unit, origin): +def write_fits(header, imagetype, data, filename, unit, origin, reorder_axis=True): """ Write a dictionary and a dataset to a FITS file Args: @@ -73,7 +73,10 @@ def write_fits(header, imagetype, data, filename, unit, origin): header['ORIGIN'] = f'Astrohack v{astrohack.__version__}: {origin}' header['DATE'] = datetime.datetime.now().strftime('%b %d %Y, %H:%M:%S') - hdu = fits.PrimaryHDU(_reorder_axes_for_fits(data)) + if reorder_axis: + hdu = fits.PrimaryHDU(_reorder_axes_for_fits(data)) + else: + hdu = fits.PrimaryHDU(data) for key in header.keys(): hdu.header.set(key, header[key]) hdu.writeto(add_prefix(filename, origin), overwrite=True) From 7b7faa5cfecceeec0ccd6eb1bac30384edc734c0 Mon Sep 17 00:00:00 2001 From: Victor de Souza magalhaes Date: Thu, 20 Feb 2025 16:42:41 -0700 Subject: [PATCH 14/53] Changed value of none clippign to the smallest value in amplitude to avoid errors when writting fits files --- src/astrohack/antenna/antenna_surface.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/astrohack/antenna/antenna_surface.py b/src/astrohack/antenna/antenna_surface.py index acd2b7fb..5f4e85c8 100644 --- a/src/astrohack/antenna/antenna_surface.py +++ b/src/astrohack/antenna/antenna_surface.py @@ -172,7 +172,7 @@ def _read_xds(self, inputxds): def _define_amp_clip(self, clip_type, clip_level): self.amplitude_noise = np.where(self.base_mask, np.nan, self.amplitude) if clip_type is None or clip_type == 'none': - clip = -np.inf + clip = np.nanmin(self.amplitude) elif clip_type == 'relative': clip = clip_level * np.nanmax(self.amplitude) elif clip_type == 'absolute': From b9672374fa677075e6913bb702e8971ae2bc3ad8 Mon Sep 17 00:00:00 2001 From: Victor de Souza magalhaes Date: Thu, 20 Feb 2025 16:43:39 -0700 Subject: [PATCH 15/53] Commented out execution code, simplified reading. --- src/astrohack/core/image_compare_tool.py | 45 ++++++++++++------------ 1 file changed, 22 insertions(+), 23 deletions(-) diff --git a/src/astrohack/core/image_compare_tool.py b/src/astrohack/core/image_compare_tool.py index 0d807bb7..7ebc2c50 100644 --- a/src/astrohack/core/image_compare_tool.py +++ b/src/astrohack/core/image_compare_tool.py @@ -22,8 +22,7 @@ class FITSImage: def __init__(self, filename, telescope_obj): self.telescope = telescope_obj - self.header, data = read_fits(filename) - self.data = data[0, 0, :, :] + self.header, self.data = read_fits(filename) self.x_axis, _, self.x_unit = get_axis_from_fits_header(self.header, 1) self.y_axis, _, self.y_unit = get_axis_from_fits_header(self.header, 2) @@ -177,27 +176,27 @@ def print_stats(self): # instatiation -first_image = image(args.first, args.noise_clip, args.blocage, - args.diameter/2, args.no_division, args.shadow_width, - args.shadow_rotation) -second_image = image(args.second, args.noise_clip, args.blocage, - args.diameter/2, args.no_division, args.shadow_width, - args.shadow_rotation) - -# Data manipulation -second_image.resample(first_image) -first_image.make_comparison(second_image) - -# Plotting -first_image.plot(args.noise_map, args.colormap, args.first_zscale) -second_image.plot(args.noise_map, args.colormap, args.second_zscale) - -if args.fits: - first_image.to_fits() - second_image.to_fits() - -if not args.quiet: - first_image.print_stats() +# first_image = image(args.first, args.noise_clip, args.blocage, +# args.diameter/2, args.no_division, args.shadow_width, +# args.shadow_rotation) +# second_image = image(args.second, args.noise_clip, args.blocage, +# args.diameter/2, args.no_division, args.shadow_width, +# args.shadow_rotation) +# +# # Data manipulation +# second_image.resample(first_image) +# first_image.make_comparison(second_image) +# +# # Plotting +# first_image.plot(args.noise_map, args.colormap, args.first_zscale) +# second_image.plot(args.noise_map, args.colormap, args.second_zscale) +# +# if args.fits: +# first_image.to_fits() +# second_image.to_fits() +# +# if not args.quiet: +# first_image.print_stats() From 2cda2d824914d7842c2f18afda3a3b9d3ed4019d Mon Sep 17 00:00:00 2001 From: Victor de Souza Magalhaes Date: Fri, 21 Feb 2025 11:06:53 -0700 Subject: [PATCH 16/53] Added a function for safe retrieval of keywords from headers/dictionaries --- src/astrohack/utils/fits.py | 20 ++++++++++++++++++-- 1 file changed, 18 insertions(+), 2 deletions(-) diff --git a/src/astrohack/utils/fits.py b/src/astrohack/utils/fits.py index e6326304..7fc073fd 100644 --- a/src/astrohack/utils/fits.py +++ b/src/astrohack/utils/fits.py @@ -10,6 +10,22 @@ from astrohack.utils.text import add_prefix +def safe_keyword_fetch(header_dict, keyword): + """ + Tries to fetch a keyword from a FITS header / dictionary + Args: + header_dict: FITS header / Dictionary + keyword: The intended keyword to fetch + + Returns: + Keyword value if prensent, None if not present. + """ + try: + return header_dict[keyword] + except KeyError: + return None + + def read_fits(filename): """ Reads a square FITS file and do sanity checks on its dimensionality @@ -51,8 +67,8 @@ def get_axis_from_fits_header(header, iaxis): val = header[f'CRVAL{iaxis}'] + inc # This makes this routine symmetrical to the put routine. axis = np.arange(n_elem) axis = val + (ref-axis)*inc - axis_unit = header[f'CUNIT{iaxis}'] - axis_type = header[f'CTYPE{iaxis}'] + axis_unit = safe_keyword_fetch(header, f'CUNIT{iaxis}') + axis_type = safe_keyword_fetch(header, f'CTYPE{iaxis}') return axis, axis_type, axis_unit From 829fd2a07c9422282d6c96f782d8d8e90578c6c2 Mon Sep 17 00:00:00 2001 From: Victor de Souza Magalhaes Date: Fri, 21 Feb 2025 11:31:56 -0700 Subject: [PATCH 17/53] Added a function to retrieve which of the axis is the stokes axis. --- src/astrohack/utils/fits.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/src/astrohack/utils/fits.py b/src/astrohack/utils/fits.py index 7fc073fd..680b02b8 100644 --- a/src/astrohack/utils/fits.py +++ b/src/astrohack/utils/fits.py @@ -10,6 +10,24 @@ from astrohack.utils.text import add_prefix +def get_stokes_axis_iaxis(header): + """ + Get which of the axis in the header is the stokes axis + Args: + header: FITS header + + Returns: + None if no stokes axis is found, iaxis if stokes axis is found + """ + naxis = header['NAXIS'] + for iaxis in range(naxis): + axis_type = safe_keyword_fetch(header, f'CTYPE{iaxis+1}') + if 'STOKES' in axis_type: + return iaxis + 1 + return None + + + def safe_keyword_fetch(header_dict, keyword): """ Tries to fetch a keyword from a FITS header / dictionary From d097464b294612b005ffa339bc6cde40ce370f8c Mon Sep 17 00:00:00 2001 From: Victor de Souza Magalhaes Date: Fri, 21 Feb 2025 11:32:28 -0700 Subject: [PATCH 18/53] Improved reading of input images --- src/astrohack/core/image_compare_tool.py | 40 +++++++++++++++++------- 1 file changed, 28 insertions(+), 12 deletions(-) diff --git a/src/astrohack/core/image_compare_tool.py b/src/astrohack/core/image_compare_tool.py index 7ebc2c50..8506f59c 100644 --- a/src/astrohack/core/image_compare_tool.py +++ b/src/astrohack/core/image_compare_tool.py @@ -7,7 +7,8 @@ from astrohack.utils.algorithms import create_aperture_mask from astrohack.visualization.plot_tools import well_positioned_colorbar from astrohack.visualization.plot_tools import close_figure, get_proper_color_map -from astrohack.utils.fits import read_fits, put_axis_in_fits_header, write_fits, get_axis_from_fits_header +from astrohack.utils.fits import read_fits, put_axis_in_fits_header, write_fits, get_axis_from_fits_header, \ + get_stokes_axis_iaxis import datetime @@ -20,15 +21,28 @@ def test_image(fits_image): class FITSImage: - def __init__(self, filename, telescope_obj): + def __init__(self, filename, telescope_obj, istokes=0, ichan=0): self.telescope = telescope_obj self.header, self.data = read_fits(filename) self.x_axis, _, self.x_unit = get_axis_from_fits_header(self.header, 1) self.y_axis, _, self.y_unit = get_axis_from_fits_header(self.header, 2) + stokes_iaxis = get_stokes_axis_iaxis(self.header) + self.unit = self.header['BUNIT'] + if len(self.data.shape) == 4: + if stokes_iaxis == 4: + self.data = self.data[istokes, ichan, ...] + else: + self.data = self.data[ichan, istokes, ...] + + elif len(self.data.shape) == 2: + pass # image is already as expected + else: + raise Exception(f'FITS image has an unsupported shape: {self.data.shape}') + if 'AIPS' in self.header['ORIGIN']: self.x_unit = 'm' self.y_unit = 'm' @@ -40,7 +54,7 @@ def __init__(self, filename, telescope_obj): self.base_mask = create_aperture_mask(self.x_axis, self.y_axis, self.telescope.inlim, self.telescope.oulim, arm_width=self.telescope.arm_shadow_width, arm_angle=self.telescope.arm_shadow_rotation) - self.rootname = '.'.join(filename.split('.')[:-1]) + self.rootname = '.'.join(filename.split('.')[:-1])+'.' self.factor = 1.0 self.residuals = None self.residuals_percent = None @@ -88,19 +102,19 @@ def image_has_same_sampling(self, ref_image): def _mask_array(self, image_array): return np.where(self.base_mask, image_array, np.nan) - def plot_results(self, destination, plot_data=False, plot_percentuals=False, plot_divided_image=False, - colormap='viridis', dpi=300, display=False): + def plot_results(self, destination, plot_residuals=True, plot_data=False, plot_percentuals=False, + plot_divided_image=False, colormap='viridis', dpi=300, display=False): extent = compute_extent(self.x_axis, self.y_axis, 0.0) cmap = get_proper_color_map(colormap) base_name = f'{destination}/{self.rootname}' - if self.residuals is None: - raise Exception("Cannot plot results as they don't exist yet.") - - self._plot_map(self._mask_array(self.residuals), 'Residuals', f'Residuals [{self.unit}]', - f'{base_name}residuals.png', cmap, extent, 'Symmetrical', dpi, display, - add_statistics=True) + if plot_residuals: + if self.residuals is None: + raise Exception("Cannot plot results as they don't exist yet.") + self._plot_map(self._mask_array(self.residuals), 'Residuals', f'Residuals [{self.unit}]', + f'{base_name}residuals.png', cmap, extent, 'symmetrical', dpi, display, + add_statistics=True) if plot_data: self._plot_map(self._mask_array(self.residuals), 'Original Data', f'Data [{self.unit}]', @@ -108,8 +122,10 @@ def plot_results(self, destination, plot_data=False, plot_percentuals=False, plo add_statistics=False) if plot_percentuals: + if self.residuals is None: + raise Exception("Cannot plot results as they don't exist yet.") self._plot_map(self._mask_array(self.residuals_percent), 'Residuals in %', f'Residuals [%]', - f'{base_name}residuals_percent.png', cmap, extent, 'Symmetrical', dpi, display, + f'{base_name}residuals_percent.png', cmap, extent, 'symmetrical', dpi, display, add_statistics=True) if plot_divided_image: From 9b4dd68a83ebb6cf0cd17f9d6812da0e3b6584ac Mon Sep 17 00:00:00 2001 From: Victor de Souza Magalhaes Date: Fri, 21 Feb 2025 11:43:41 -0700 Subject: [PATCH 19/53] Added control over pixel_offset being applied --- src/astrohack/utils/fits.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/astrohack/utils/fits.py b/src/astrohack/utils/fits.py index 680b02b8..d72ef301 100644 --- a/src/astrohack/utils/fits.py +++ b/src/astrohack/utils/fits.py @@ -69,12 +69,13 @@ def read_fits(filename): return head, data -def get_axis_from_fits_header(header, iaxis): +def get_axis_from_fits_header(header, iaxis, pixel_offset=True): """ Pull axis information from FITS file and store it in a numpy array, ignores rotation in axes. Args: header: FITS header iaxis: Which axis is to be fetched from the header. + pixel_offset: apply one pixel offset Returns: numpy array representation of axis, axis type and axis unit @@ -82,7 +83,10 @@ def get_axis_from_fits_header(header, iaxis): n_elem = header[f'NAXIS{iaxis}'] ref = header[f'CRPIX{iaxis}'] inc = header[f'CDELT{iaxis}'] - val = header[f'CRVAL{iaxis}'] + inc # This makes this routine symmetrical to the put routine. + if pixel_offset: + val = header[f'CRVAL{iaxis}'] + inc # This makes this routine symmetrical to the put routine. + else: + val = header[f'CRVAL{iaxis}'] axis = np.arange(n_elem) axis = val + (ref-axis)*inc axis_unit = safe_keyword_fetch(header, f'CUNIT{iaxis}') From 6d29d73f557e627ec1087b596d2f0e555f7b7f3a Mon Sep 17 00:00:00 2001 From: Victor de Souza Magalhaes Date: Fri, 21 Feb 2025 11:44:09 -0700 Subject: [PATCH 20/53] Fixed data ordering and adapted axis reading --- src/astrohack/core/image_compare_tool.py | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/src/astrohack/core/image_compare_tool.py b/src/astrohack/core/image_compare_tool.py index 8506f59c..64d305f3 100644 --- a/src/astrohack/core/image_compare_tool.py +++ b/src/astrohack/core/image_compare_tool.py @@ -25,9 +25,6 @@ def __init__(self, filename, telescope_obj, istokes=0, ichan=0): self.telescope = telescope_obj self.header, self.data = read_fits(filename) - self.x_axis, _, self.x_unit = get_axis_from_fits_header(self.header, 1) - self.y_axis, _, self.y_unit = get_axis_from_fits_header(self.header, 2) - stokes_iaxis = get_stokes_axis_iaxis(self.header) self.unit = self.header['BUNIT'] @@ -44,10 +41,16 @@ def __init__(self, filename, telescope_obj, istokes=0, ichan=0): raise Exception(f'FITS image has an unsupported shape: {self.data.shape}') if 'AIPS' in self.header['ORIGIN']: + self.x_axis, _, self.x_unit = get_axis_from_fits_header(self.header, 1, pixel_offset=False) + self.y_axis, _, self.y_unit = get_axis_from_fits_header(self.header, 2, pixel_offset=False) + self.x_unit = 'm' self.y_unit = 'm' elif 'Astrohack' in self.header['ORIGIN']: - self.data = self.data.T + self.x_axis, _, self.x_unit = get_axis_from_fits_header(self.header, 1) + self.y_axis, _, self.y_unit = get_axis_from_fits_header(self.header, 2) + + self.data = np.fliplr(self.data) else: raise Exception(f'Unrecognized origin:\n{self.header["origin"]}') @@ -102,8 +105,8 @@ def image_has_same_sampling(self, ref_image): def _mask_array(self, image_array): return np.where(self.base_mask, image_array, np.nan) - def plot_results(self, destination, plot_residuals=True, plot_data=False, plot_percentuals=False, - plot_divided_image=False, colormap='viridis', dpi=300, display=False): + def plot_images(self, destination, plot_residuals=True, plot_data=False, plot_percentuals=False, + plot_divided_image=False, colormap='viridis', dpi=300, display=False): extent = compute_extent(self.x_axis, self.y_axis, 0.0) cmap = get_proper_color_map(colormap) @@ -117,7 +120,7 @@ def plot_results(self, destination, plot_residuals=True, plot_data=False, plot_p add_statistics=True) if plot_data: - self._plot_map(self._mask_array(self.residuals), 'Original Data', f'Data [{self.unit}]', + self._plot_map(self._mask_array(self.data), 'Original Data', f'Data [{self.unit}]', f'{base_name}data.png', cmap, extent, [None, None], dpi, display, add_statistics=False) From fd2fa83a8c44c0f0b8144c81d7e5572086f634a1 Mon Sep 17 00:00:00 2001 From: Victor de Souza Magalhaes Date: Fri, 21 Feb 2025 15:25:05 -0700 Subject: [PATCH 21/53] Added xds export --- src/astrohack/core/image_compare_tool.py | 25 ++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/src/astrohack/core/image_compare_tool.py b/src/astrohack/core/image_compare_tool.py index 64d305f3..96f1bd87 100644 --- a/src/astrohack/core/image_compare_tool.py +++ b/src/astrohack/core/image_compare_tool.py @@ -1,6 +1,7 @@ import numpy as np from scipy.interpolate import griddata from matplotlib import pyplot as plt +import xarray as xr from astrohack import compute_extent, data_statistics, statistics_to_text from astrohack.utils import are_axes_equal @@ -160,6 +161,30 @@ def _plot_map(self, data, title, zlabel, filename, cmap, extent, zscale, dpi, di ax.set_title(statistics_to_text(data_stats)) close_figure(fig, title, filename, dpi, display) + def export_as_xds(self): + xds = xr.Dataset() + obj_dict = vars(self) + + coords = {'x': self.x_axis, 'y': self.y_axis} + for key, value in obj_dict.items(): + failed = False + if isinstance(value, np.ndarray): + if len(value.shape) == 2: + xds[key] = xr.DataArray(value, dims=['x', 'y']) + elif len(value.shape) == 1: + pass # Axes + else: + failed = True + else: + xds.attrs[key] = value + + if failed: + raise Exception(f"Don't know what to do with: {key}") + xds.assign_coords(coords) + return xds + + + # def to_fits(self): # fits = '.fits' # header = self.header.copy() From 373270d8c556e9265a692cad47e8730703a1ee6f Mon Sep 17 00:00:00 2001 From: Victor de Souza Magalhaes Date: Fri, 21 Feb 2025 15:55:47 -0700 Subject: [PATCH 22/53] Added the option to read a FITS header as a python dictionary --- src/astrohack/utils/fits.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/src/astrohack/utils/fits.py b/src/astrohack/utils/fits.py index d72ef301..b6ff7c0e 100644 --- a/src/astrohack/utils/fits.py +++ b/src/astrohack/utils/fits.py @@ -44,7 +44,7 @@ def safe_keyword_fetch(header_dict, keyword): return None -def read_fits(filename): +def read_fits(filename, header_as_dict=True): """ Reads a square FITS file and do sanity checks on its dimensionality Args: @@ -66,7 +66,14 @@ def read_fits(filename): raise Exception(filename + " is not bi-dimensional") if head["NAXIS1"] != head["NAXIS2"]: raise Exception(filename + " does not have the same amount of pixels in the x and y axes") - return head, data + + if header_as_dict: + header_dict = {} + for key, value in head.items(): + header_dict[key] = value + return header_dict, data + else: + return head, data def get_axis_from_fits_header(header, iaxis, pixel_offset=True): From b1ad5ca7cdaaf314b17ad8d859f2ab0f68b32493 Mon Sep 17 00:00:00 2001 From: Victor de Souza Magalhaes Date: Fri, 21 Feb 2025 16:22:23 -0700 Subject: [PATCH 23/53] Factorised mask creation, added read and write to Zarr and a print method --- src/astrohack/core/image_compare_tool.py | 111 ++++++++++++++++------- 1 file changed, 79 insertions(+), 32 deletions(-) diff --git a/src/astrohack/core/image_compare_tool.py b/src/astrohack/core/image_compare_tool.py index 96f1bd87..7109129c 100644 --- a/src/astrohack/core/image_compare_tool.py +++ b/src/astrohack/core/image_compare_tool.py @@ -3,14 +3,12 @@ from matplotlib import pyplot as plt import xarray as xr -from astrohack import compute_extent, data_statistics, statistics_to_text -from astrohack.utils import are_axes_equal -from astrohack.utils.algorithms import create_aperture_mask -from astrohack.visualization.plot_tools import well_positioned_colorbar +from astrohack.antenna.telescope import Telescope +from astrohack.utils.text import statistics_to_text +from astrohack.utils.algorithms import create_aperture_mask, data_statistics, are_axes_equal +from astrohack.visualization.plot_tools import well_positioned_colorbar, compute_extent from astrohack.visualization.plot_tools import close_figure, get_proper_color_map -from astrohack.utils.fits import read_fits, put_axis_in_fits_header, write_fits, get_axis_from_fits_header, \ - get_stokes_axis_iaxis -import datetime +from astrohack.utils.fits import read_fits, get_axis_from_fits_header, get_stokes_axis_iaxis def test_image(fits_image): @@ -22,10 +20,38 @@ def test_image(fits_image): class FITSImage: - def __init__(self, filename, telescope_obj, istokes=0, ichan=0): - self.telescope = telescope_obj - self.header, self.data = read_fits(filename) + def __init__(self, filename: str, telescope_name: str): + # Initialization from parameters + self.filename = filename + self.telescope_name = telescope_name + self.rootname = '.'.join(filename.split('.')[:-1])+'.' + + # Blank slate initialization + self.header = None + self.data = None + self.factor = 1.0 + self.residuals = None + self.residuals_percent = None + self.divided_image = None + self.reference_name = None + self.resampled = False + self.x_axis = None + self.y_axis = None + self.x_unit = None + self.y_unit = None + self.unit = None + self.fits_name = None + + if '.FITS' in filename.upper(): + self._init_as_fits(0, 0) + elif '.zarr' in filename: + self._init_as_xds() + else: + raise Exception(f"Don't know how to read {filename}") + def _init_as_fits(self, istokes, ichan): + self.header, self.data = read_fits(self.filename, header_as_dict=True) + self.fits_name = self.filename stokes_iaxis = get_stokes_axis_iaxis(self.header) self.unit = self.header['BUNIT'] @@ -44,26 +70,32 @@ def __init__(self, filename, telescope_obj, istokes=0, ichan=0): if 'AIPS' in self.header['ORIGIN']: self.x_axis, _, self.x_unit = get_axis_from_fits_header(self.header, 1, pixel_offset=False) self.y_axis, _, self.y_unit = get_axis_from_fits_header(self.header, 2, pixel_offset=False) - self.x_unit = 'm' self.y_unit = 'm' elif 'Astrohack' in self.header['ORIGIN']: self.x_axis, _, self.x_unit = get_axis_from_fits_header(self.header, 1) self.y_axis, _, self.y_unit = get_axis_from_fits_header(self.header, 2) - self.data = np.fliplr(self.data) else: raise Exception(f'Unrecognized origin:\n{self.header["origin"]}') + self._create_base_mask() - self.base_mask = create_aperture_mask(self.x_axis, self.y_axis, self.telescope.inlim, self.telescope.oulim, - arm_width=self.telescope.arm_shadow_width, - arm_angle=self.telescope.arm_shadow_rotation) - self.rootname = '.'.join(filename.split('.')[:-1])+'.' - self.factor = 1.0 - self.residuals = None - self.residuals_percent = None - self.divided_image = None - self.reference_name = None + def _init_as_xds(self): + xds = xr.open_zarr(self.filename) + for key in xds.attrs: + setattr(self, key, xds.attrs[key]) + + self.x_axis = xds.x.values + self.y_axis = xds.y.values + + for key, value in xds.items(): + setattr(self, key, xds[key].values) + + def _create_base_mask(self): + telescope_obj = Telescope(self.telescope_name) + self.base_mask = create_aperture_mask(self.x_axis, self.y_axis, telescope_obj.inlim, telescope_obj.oulim, + arm_width=telescope_obj.arm_shadow_width, + arm_angle=telescope_obj.arm_shadow_rotation) def resample(self, ref_image): test_image(ref_image) @@ -76,9 +108,8 @@ def resample(self, ref_image): self.x_axis = ref_image.x_axis self.y_axis = ref_image.y_axis self.data = resamp.reshape(size) - self.base_mask = create_aperture_mask(self.x_axis, self.y_axis, self.telescope.inlim, self.telescope.oulim, - arm_width=self.telescope.arm_shadow_width, - arm_angle=self.telescope.arm_shadow_rotation) + self._create_base_mask() + self.resampled = True def compare_difference(self, ref_image): test_image(ref_image) @@ -180,10 +211,26 @@ def export_as_xds(self): if failed: raise Exception(f"Don't know what to do with: {key}") - xds.assign_coords(coords) + + xds = xds.assign_coords(coords) return xds + def to_zarr(self, zarr_filename): + xds = self.export_as_xds() + xds.to_zarr(zarr_filename, mode="w", compute=True, consolidated=True) + def __repr__(self): + obj_dict = vars(self) + outstr = '' + for key, value in obj_dict.items(): + if isinstance(value, np.ndarray): + outstr += f'{key:17s} -> {value.shape}' + elif isinstance(value, dict): + outstr += f'{key:17s} -> dict()' + else: + outstr += f'{key:17s} = {value}' + outstr += '\n' + return outstr # def to_fits(self): # fits = '.fits' @@ -210,13 +257,13 @@ def export_as_xds(self): # create_fits(header, self.noise, f'{filename}.noise{fits}') - def print_stats(self): - print(80*'*') - print() - print(f'Mean scaling factor = {self.factor:.3}') - print(f'Mean Residual = {self.res_mean:.3}%') - print(f'Residuals RMS = {self.res_rms:.3}%') - print() + # def print_stats(self): + # print(80*'*') + # print() + # print(f'Mean scaling factor = {self.factor:.3}') + # print(f'Mean Residual = {self.res_mean:.3}%') + # print(f'Residuals RMS = {self.res_rms:.3}%') + # print() # instatiation From 22ba437091292093cc1785e06ebd0478bc854e09 Mon Sep 17 00:00:00 2001 From: Victor de Souza magalhaes Date: Tue, 25 Feb 2025 08:55:41 -0700 Subject: [PATCH 24/53] Input file name is preserved when reading a .zarr container. --- src/astrohack/core/image_compare_tool.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/astrohack/core/image_compare_tool.py b/src/astrohack/core/image_compare_tool.py index 7109129c..26661bac 100644 --- a/src/astrohack/core/image_compare_tool.py +++ b/src/astrohack/core/image_compare_tool.py @@ -81,6 +81,7 @@ def _init_as_fits(self, istokes, ichan): self._create_base_mask() def _init_as_xds(self): + filename = self.filename xds = xr.open_zarr(self.filename) for key in xds.attrs: setattr(self, key, xds.attrs[key]) @@ -91,6 +92,8 @@ def _init_as_xds(self): for key, value in xds.items(): setattr(self, key, xds[key].values) + self.filename = filename + def _create_base_mask(self): telescope_obj = Telescope(self.telescope_name) self.base_mask = create_aperture_mask(self.x_axis, self.y_axis, telescope_obj.inlim, telescope_obj.oulim, From 2322c58a8aa35b7f7339665ac4d699e67498cb8e Mon Sep 17 00:00:00 2001 From: Victor de Souza magalhaes Date: Tue, 25 Feb 2025 09:00:50 -0700 Subject: [PATCH 25/53] Added storage of original data and axes description --- src/astrohack/core/image_compare_tool.py | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/src/astrohack/core/image_compare_tool.py b/src/astrohack/core/image_compare_tool.py index 26661bac..b4c8de14 100644 --- a/src/astrohack/core/image_compare_tool.py +++ b/src/astrohack/core/image_compare_tool.py @@ -37,10 +37,13 @@ def __init__(self, filename: str, telescope_name: str): self.resampled = False self.x_axis = None self.y_axis = None + self.original_x_axis = None + self.original_y_axis = None self.x_unit = None self.y_unit = None self.unit = None self.fits_name = None + self.original_data = None if '.FITS' in filename.upper(): self._init_as_fits(0, 0) @@ -51,6 +54,7 @@ def __init__(self, filename: str, telescope_name: str): def _init_as_fits(self, istokes, ichan): self.header, self.data = read_fits(self.filename, header_as_dict=True) + self.original_data = np.copy(self.data) self.fits_name = self.filename stokes_iaxis = get_stokes_axis_iaxis(self.header) @@ -63,7 +67,7 @@ def _init_as_fits(self, istokes, ichan): self.data = self.data[ichan, istokes, ...] elif len(self.data.shape) == 2: - pass # image is already as expected + pass # image is already as expected else: raise Exception(f'FITS image has an unsupported shape: {self.data.shape}') @@ -79,6 +83,8 @@ def _init_as_fits(self, istokes, ichan): else: raise Exception(f'Unrecognized origin:\n{self.header["origin"]}') self._create_base_mask() + self.original_x_axis = np.copy(self.x_axis) + self.original_y_axis = np.copy(self.y_axis) def _init_as_xds(self): filename = self.filename @@ -88,6 +94,8 @@ def _init_as_xds(self): self.x_axis = xds.x.values self.y_axis = xds.y.values + self.original_x_axis = xds.original_x.values + self.original_y_axis = xds.original_y.values for key, value in xds.items(): setattr(self, key, xds[key].values) @@ -199,12 +207,16 @@ def export_as_xds(self): xds = xr.Dataset() obj_dict = vars(self) - coords = {'x': self.x_axis, 'y': self.y_axis} + coords = {'x': self.x_axis, 'y': self.y_axis, + 'original_x': self.original_x_axis, 'original_y': self.original_y_axis} for key, value in obj_dict.items(): failed = False if isinstance(value, np.ndarray): if len(value.shape) == 2: - xds[key] = xr.DataArray(value, dims=['x', 'y']) + if 'original' in key: + xds[key] = xr.DataArray(value, dims=['original_x', 'original_y']) + else: + xds[key] = xr.DataArray(value, dims=['x', 'y']) elif len(value.shape) == 1: pass # Axes else: From 713ad0664cb5072b86b2604997ee87abdc444434 Mon Sep 17 00:00:00 2001 From: Victor de Souza Magalhaes Date: Tue, 25 Feb 2025 15:39:58 -0700 Subject: [PATCH 26/53] Made FITS write more general --- src/astrohack/utils/fits.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/src/astrohack/utils/fits.py b/src/astrohack/utils/fits.py index b6ff7c0e..f0d8d5b3 100644 --- a/src/astrohack/utils/fits.py +++ b/src/astrohack/utils/fits.py @@ -101,7 +101,7 @@ def get_axis_from_fits_header(header, iaxis, pixel_offset=True): return axis, axis_type, axis_unit -def write_fits(header, imagetype, data, filename, unit, origin, reorder_axis=True): +def write_fits(header, imagetype, data, filename, unit, origin=None, reorder_axis=True): """ Write a dictionary and a dataset to a FITS file Args: @@ -111,6 +111,7 @@ def write_fits(header, imagetype, data, filename, unit, origin, reorder_axis=Tru filename: The name of the output file unit: to be set to bunit origin: Which astrohack mds has created the FITS being written + reorder_axis: Reorder data axes so that they are compatible with regular FITS ordering """ header['BUNIT'] = unit @@ -118,13 +119,20 @@ def write_fits(header, imagetype, data, filename, unit, origin, reorder_axis=Tru header['ORIGIN'] = f'Astrohack v{astrohack.__version__}: {origin}' header['DATE'] = datetime.datetime.now().strftime('%b %d %Y, %H:%M:%S') + if origin is None: + header['ORIGIN'] = f'Astrohack v{astrohack.__version__}' + outfile = filename + else: + header['ORIGIN'] = f'Astrohack v{astrohack.__version__}: {origin}' + outfile = add_prefix(filename, origin) + if reorder_axis: hdu = fits.PrimaryHDU(_reorder_axes_for_fits(data)) else: hdu = fits.PrimaryHDU(data) for key in header.keys(): hdu.header.set(key, header[key]) - hdu.writeto(add_prefix(filename, origin), overwrite=True) + hdu.writeto(outfile, overwrite=True) return From 8caf278080f30ced76d562dbdc4f06c2d0b3ee1a Mon Sep 17 00:00:00 2001 From: Victor de Souza Magalhaes Date: Tue, 25 Feb 2025 15:40:16 -0700 Subject: [PATCH 27/53] Added method to export to FITS --- src/astrohack/core/image_compare_tool.py | 50 ++++++++++++------------ 1 file changed, 26 insertions(+), 24 deletions(-) diff --git a/src/astrohack/core/image_compare_tool.py b/src/astrohack/core/image_compare_tool.py index b4c8de14..fc00bd65 100644 --- a/src/astrohack/core/image_compare_tool.py +++ b/src/astrohack/core/image_compare_tool.py @@ -2,13 +2,15 @@ from scipy.interpolate import griddata from matplotlib import pyplot as plt import xarray as xr +import pathlib from astrohack.antenna.telescope import Telescope from astrohack.utils.text import statistics_to_text from astrohack.utils.algorithms import create_aperture_mask, data_statistics, are_axes_equal from astrohack.visualization.plot_tools import well_positioned_colorbar, compute_extent from astrohack.visualization.plot_tools import close_figure, get_proper_color_map -from astrohack.utils.fits import read_fits, get_axis_from_fits_header, get_stokes_axis_iaxis +from astrohack.utils.fits import read_fits, get_axis_from_fits_header, get_stokes_axis_iaxis, put_axis_in_fits_header, \ + write_fits def test_image(fits_image): @@ -247,30 +249,30 @@ def __repr__(self): outstr += '\n' return outstr - # def to_fits(self): - # fits = '.fits' - # header = self.header.copy() - # put_axis_in_header(self.x_axis, self.x_unit, 1, header) - # put_axis_in_header(self.y_axis, self.y_unit, 2, header) - # - # if self.resampled: - # filename = f'comp_{self.rootname}.resampled' - # else: - # filename = f'comp_{self.rootname}' - # - # create_fits(header, self.masked, f'{filename}.masked{fits}') - # create_fits(header, self.mask, f'{filename}.mask{fits}') - # - # if self.division is not None: - # create_fits(header, self.division, f'{filename}.division{fits}') - # - # if self.residuals is not None: - # create_fits(header, self.residuals, f'{filename}.residual{fits}') - # - # - # if self.noise is not None: - # create_fits(header, self.noise, f'{filename}.noise{fits}') + def export_to_fits(self, destination): + pathlib.Path(destination).mkdir(exist_ok=True) + ext_fits = '.fits' + out_header = self.header.copy() + put_axis_in_fits_header(out_header, self.x_axis, 1, '', self.x_unit) + put_axis_in_fits_header(out_header, self.y_axis, 2, '', self.y_unit) + + obj_dict = vars(self) + for key, value in obj_dict.items(): + if isinstance(value, np.ndarray): + if len(value.shape) == 2: + if 'original' in key: + pass + else: + if key == 'base_mask' or key == 'divided_image': + unit = '' + + elif key == 'residuals_percent': + unit = '%' + else: + unit = self.unit + filename = f'{destination}/{self.rootname}{key}{ext_fits}' + write_fits(out_header, key, np.fliplr(value.astype(float)), filename, unit, reorder_axis=False) # def print_stats(self): # print(80*'*') From 42eb85f836310e5a9f1c4d30b97adecd67a1238d Mon Sep 17 00:00:00 2001 From: Victor de Souza Magalhaes Date: Tue, 25 Feb 2025 15:47:00 -0700 Subject: [PATCH 28/53] Cleanup --- src/astrohack/core/image_compare_tool.py | 41 ++---------------------- 1 file changed, 2 insertions(+), 39 deletions(-) diff --git a/src/astrohack/core/image_compare_tool.py b/src/astrohack/core/image_compare_tool.py index fc00bd65..7e60112c 100644 --- a/src/astrohack/core/image_compare_tool.py +++ b/src/astrohack/core/image_compare_tool.py @@ -48,13 +48,13 @@ def __init__(self, filename: str, telescope_name: str): self.original_data = None if '.FITS' in filename.upper(): - self._init_as_fits(0, 0) + self._init_as_fits() elif '.zarr' in filename: self._init_as_xds() else: raise Exception(f"Don't know how to read {filename}") - def _init_as_fits(self, istokes, ichan): + def _init_as_fits(self, istokes=0, ichan=0): self.header, self.data = read_fits(self.filename, header_as_dict=True) self.original_data = np.copy(self.data) self.fits_name = self.filename @@ -273,40 +273,3 @@ def export_to_fits(self, destination): unit = self.unit filename = f'{destination}/{self.rootname}{key}{ext_fits}' write_fits(out_header, key, np.fliplr(value.astype(float)), filename, unit, reorder_axis=False) - - # def print_stats(self): - # print(80*'*') - # print() - # print(f'Mean scaling factor = {self.factor:.3}') - # print(f'Mean Residual = {self.res_mean:.3}%') - # print(f'Residuals RMS = {self.res_rms:.3}%') - # print() - - -# instatiation -# first_image = image(args.first, args.noise_clip, args.blocage, -# args.diameter/2, args.no_division, args.shadow_width, -# args.shadow_rotation) -# second_image = image(args.second, args.noise_clip, args.blocage, -# args.diameter/2, args.no_division, args.shadow_width, -# args.shadow_rotation) -# -# # Data manipulation -# second_image.resample(first_image) -# first_image.make_comparison(second_image) -# -# # Plotting -# first_image.plot(args.noise_map, args.colormap, args.first_zscale) -# second_image.plot(args.noise_map, args.colormap, args.second_zscale) -# -# if args.fits: -# first_image.to_fits() -# second_image.to_fits() -# -# if not args.quiet: -# first_image.print_stats() - - - - - From 16a04bbaf2e2b4b8e7ce51e450ddfd4d1e67c2fe Mon Sep 17 00:00:00 2001 From: Victor de Souza magalhaes Date: Wed, 26 Feb 2025 14:15:07 -0700 Subject: [PATCH 29/53] Minor documentation. --- src/astrohack/core/image_compare_tool.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/astrohack/core/image_compare_tool.py b/src/astrohack/core/image_compare_tool.py index 7e60112c..cea9d48a 100644 --- a/src/astrohack/core/image_compare_tool.py +++ b/src/astrohack/core/image_compare_tool.py @@ -23,6 +23,13 @@ def test_image(fits_image): class FITSImage: def __init__(self, filename: str, telescope_name: str): + """ + Initializes the FITSImage object from a file on disk + Args: + filename: Name of the file on disk, may be .FITS file or a .zarr xds with a disk representation of a \ + FITSImage object + telescope_name: Name of the telescope used on the images so that masking can be properly applied. + """ # Initialization from parameters self.filename = filename self.telescope_name = telescope_name From a4aebb2527ddfac475b850cae107cf9aaffd6620 Mon Sep 17 00:00:00 2001 From: Victor de Souza magalhaes Date: Wed, 26 Feb 2025 14:27:57 -0700 Subject: [PATCH 30/53] Added a method to create a scatter plot against the reference image --- src/astrohack/core/image_compare_tool.py | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/src/astrohack/core/image_compare_tool.py b/src/astrohack/core/image_compare_tool.py index cea9d48a..43abcda2 100644 --- a/src/astrohack/core/image_compare_tool.py +++ b/src/astrohack/core/image_compare_tool.py @@ -8,7 +8,7 @@ from astrohack.utils.text import statistics_to_text from astrohack.utils.algorithms import create_aperture_mask, data_statistics, are_axes_equal from astrohack.visualization.plot_tools import well_positioned_colorbar, compute_extent -from astrohack.visualization.plot_tools import close_figure, get_proper_color_map +from astrohack.visualization.plot_tools import close_figure, get_proper_color_map, scatter_plot from astrohack.utils.fits import read_fits, get_axis_from_fits_header, get_stokes_axis_iaxis, put_axis_in_fits_header, \ write_fits @@ -280,3 +280,16 @@ def export_to_fits(self, destination): unit = self.unit filename = f'{destination}/{self.rootname}{key}{ext_fits}' write_fits(out_header, key, np.fliplr(value.astype(float)), filename, unit, reorder_axis=False) + + def scatter_plot(self, ref_image, dpi=300, display=False): + test_image(ref_image) + if not self.image_has_same_sampling(ref_image): + self.resample(ref_image) + + fig, ax = plt.subplots(1, 1, figsize=[10, 8]) + ydata = self.data[self.base_mask] + xdata = ref_image.data[self.base_mask] + scatter_plot(ax, xdata, f'Reference image {ref_image.filename} [{ref_image.unit}]', + ydata, f'{self.filename} [{self.unit}]') + close_figure(fig, 'Scatter plot against reference image', f'{self.rootname}scatter.png', dpi, display) + From 9dc627b09984669250dc122fc857985917b216b2 Mon Sep 17 00:00:00 2001 From: Victor de Souza magalhaes Date: Wed, 26 Feb 2025 14:52:17 -0700 Subject: [PATCH 31/53] Added capability to plot a linear regression together with a scatter plot --- src/astrohack/visualization/plot_tools.py | 27 ++++++++++++++++++++--- 1 file changed, 24 insertions(+), 3 deletions(-) diff --git a/src/astrohack/visualization/plot_tools.py b/src/astrohack/visualization/plot_tools.py index 618132ea..cf361a93 100644 --- a/src/astrohack/visualization/plot_tools.py +++ b/src/astrohack/visualization/plot_tools.py @@ -1,5 +1,6 @@ import matplotlib.image import numpy as np +from scipy.stats import linregress from matplotlib import pyplot as plt from matplotlib.patches import Rectangle @@ -218,7 +219,10 @@ def scatter_plot( residuals_marker='+', residuals_color='black', residuals_linestyle='', - residuals_label='residuals' + residuals_label='residuals', + add_regression=False, + regression_linestyle='-', + regression_color='black' ): """ Do scatter simple scatter plots of data to a plotting axis @@ -245,12 +249,19 @@ def scatter_plot( model_color: Color of the model marker model_linestyle: Line style for connecting model points model_label: Label for model points + plot_residuals: Add a residuals subplot at the bottom when a model is provided + residuals_marker: Marker for residuals + residuals_color: Color for residual markers + residuals_linestyle: Line style for residuals + residuals_label: Label for residuals + add_regression: Add a linear regression between X and y data + regression_linestyle: Line style for the regression plot + regression_color: Color for the regression plot """ ax.plot(xdata, ydata, ls=data_linestyle, marker=data_marker, color=data_color, label=data_label) ax.set_xlabel(xlabel) ax.set_ylabel(ylabel) - if title is not None: - ax.set_title(title) + if xlim is not None: ax.set_xlim(xlim) if ylim is not None: @@ -274,6 +285,13 @@ def scatter_plot( rotation=20 ) + if add_regression: + slope, intercept, _, _, _ = linregress(xdata, ydata) + regression_label = f'y = {slope:.4f}*x + {intercept:.4f}' + yregress = slope*xdata + intercept + ax.plot(xdata, yregress, ls=regression_linestyle, color=regression_color, label=regression_label, lw=2) + ax.legend() + if model is not None: ax.plot(xdata, model, ls=model_linestyle, marker=model_marker, color=model_color, label=model_label) ax.legend() @@ -297,4 +315,7 @@ def scatter_plot( ax_res.axhline(0, color=hv_color, ls=hv_linestyle) ax_res.set_ylabel('Residuals') + if title is not None: + ax.set_title(title) + return From 9585fba032bc2c94078d523518435808ee756a08 Mon Sep 17 00:00:00 2001 From: Victor de Souza magalhaes Date: Wed, 26 Feb 2025 14:59:52 -0700 Subject: [PATCH 32/53] Small change in formating --- src/astrohack/utils/text.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/astrohack/utils/text.py b/src/astrohack/utils/text.py index 8522386a..fecbeea9 100644 --- a/src/astrohack/utils/text.py +++ b/src/astrohack/utils/text.py @@ -574,7 +574,7 @@ def significant_figures_round(x, digits): def statistics_to_text(data_statistics, keys=None): if keys is None: - outstr = (f'min={data_statistics["min"]:.2e}, max={data_statistics["max"]:.2f}, ' + outstr = (f'min={data_statistics["min"]:.2f}, max={data_statistics["max"]:.2f}, ' f'mean={data_statistics["mean"]:.2f}, med={data_statistics["median"]:.2f}, ' f'rms={data_statistics["rms"]:.2f}') else: From 4bca3fd60d5e0c992663a724868e4272729991aa Mon Sep 17 00:00:00 2001 From: Victor de Souza magalhaes Date: Wed, 26 Feb 2025 15:04:04 -0700 Subject: [PATCH 33/53] Plotting improvements --- src/astrohack/core/image_compare_tool.py | 43 +++++++++++++----------- 1 file changed, 23 insertions(+), 20 deletions(-) diff --git a/src/astrohack/core/image_compare_tool.py b/src/astrohack/core/image_compare_tool.py index 43abcda2..a82bf3e8 100644 --- a/src/astrohack/core/image_compare_tool.py +++ b/src/astrohack/core/image_compare_tool.py @@ -138,7 +138,7 @@ def compare_difference(self, ref_image): self.residuals = ref_image.data - (self.data * self.factor) self.residuals_percent = 100 * self.residuals/ref_image.data - self.reference_name = ref_image.rootname + self.reference_name = ref_image.filename def compare_scaled_difference(self, ref_image, rejection=10): test_image(ref_image) @@ -157,19 +157,18 @@ def image_has_same_sampling(self, ref_image): def _mask_array(self, image_array): return np.where(self.base_mask, image_array, np.nan) - def plot_images(self, destination, plot_residuals=True, plot_data=False, plot_percentuals=False, + def plot_images(self, destination, plot_data=False, plot_percentuals=False, plot_divided_image=False, colormap='viridis', dpi=300, display=False): extent = compute_extent(self.x_axis, self.y_axis, 0.0) cmap = get_proper_color_map(colormap) base_name = f'{destination}/{self.rootname}' - if plot_residuals: - if self.residuals is None: - raise Exception("Cannot plot results as they don't exist yet.") - self._plot_map(self._mask_array(self.residuals), 'Residuals', f'Residuals [{self.unit}]', - f'{base_name}residuals.png', cmap, extent, 'symmetrical', dpi, display, - add_statistics=True) + if self.residuals is None: + raise Exception("Cannot plot results as they don't exist yet.") + self._plot_map(self._mask_array(self.residuals), f'Residuals, ref={self.reference_name}', + f'Residuals [{self.unit}]', f'{base_name}residuals.png', cmap, extent, + 'symmetrical', dpi, display, add_statistics=True) if plot_data: self._plot_map(self._mask_array(self.data), 'Original Data', f'Data [{self.unit}]', @@ -179,16 +178,16 @@ def plot_images(self, destination, plot_residuals=True, plot_data=False, plot_pe if plot_percentuals: if self.residuals is None: raise Exception("Cannot plot results as they don't exist yet.") - self._plot_map(self._mask_array(self.residuals_percent), 'Residuals in %', f'Residuals [%]', - f'{base_name}residuals_percent.png', cmap, extent, 'symmetrical', dpi, display, - add_statistics=True) + self._plot_map(self._mask_array(self.residuals_percent), f'Residuals in %, ref={self.reference_name}', + f'Residuals [%]', f'{base_name}residuals_percent.png', cmap, extent, + 'symmetrical', dpi, display, add_statistics=True) if plot_divided_image: if self.divided_image is None: raise Exception("Cannot plot a divided image that does not exist.") - self._plot_map(self._mask_array(self.divided_image), 'Divided image', f'Division [ ]', - f'{base_name}divided.png', cmap, extent, [None, None], dpi, display, - add_statistics=True) + self._plot_map(self._mask_array(self.divided_image), f'Divided image, ref={self.reference_name}', + f'Division [ ]', f'{base_name}divided.png', cmap, extent, [None, None], + dpi, display, add_statistics=True) def _plot_map(self, data, title, zlabel, filename, cmap, extent, zscale, dpi, display, add_statistics=False): fig, ax = plt.subplots(1, 1, figsize=[10, 8]) @@ -281,15 +280,19 @@ def export_to_fits(self, destination): filename = f'{destination}/{self.rootname}{key}{ext_fits}' write_fits(out_header, key, np.fliplr(value.astype(float)), filename, unit, reorder_axis=False) - def scatter_plot(self, ref_image, dpi=300, display=False): + def scatter_plot(self, destination, ref_image, dpi=300, display=False): test_image(ref_image) if not self.image_has_same_sampling(ref_image): self.resample(ref_image) fig, ax = plt.subplots(1, 1, figsize=[10, 8]) - ydata = self.data[self.base_mask] - xdata = ref_image.data[self.base_mask] - scatter_plot(ax, xdata, f'Reference image {ref_image.filename} [{ref_image.unit}]', - ydata, f'{self.filename} [{self.unit}]') - close_figure(fig, 'Scatter plot against reference image', f'{self.rootname}scatter.png', dpi, display) + scatter_mask = np.isfinite(ref_image.data) + scatter_mask = np.where(np.isfinite(self.data), scatter_mask, False) + ydata = self.data[scatter_mask] + xdata = ref_image.data[scatter_mask] + + scatter_plot(ax, xdata, f'Reference image {ref_image.filename} [{ref_image.unit}]', + ydata, f'{self.filename} [{self.unit}]', add_regression=True) + close_figure(fig, 'Scatter plot against reference image', f'{destination}/{self.rootname}scatter.png', + dpi, display) From 27966f34d295651c9939b327de702a84a64ba289 Mon Sep 17 00:00:00 2001 From: Victor de Souza magalhaes Date: Wed, 26 Feb 2025 15:04:39 -0700 Subject: [PATCH 34/53] Minor issues --- src/astrohack/core/image_compare_tool.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/astrohack/core/image_compare_tool.py b/src/astrohack/core/image_compare_tool.py index a82bf3e8..6859122d 100644 --- a/src/astrohack/core/image_compare_tool.py +++ b/src/astrohack/core/image_compare_tool.py @@ -107,7 +107,7 @@ def _init_as_xds(self): self.original_y_axis = xds.original_y.values for key, value in xds.items(): - setattr(self, key, xds[key].values) + setattr(self, str(key), xds[key].values) self.filename = filename @@ -226,7 +226,7 @@ def export_as_xds(self): else: xds[key] = xr.DataArray(value, dims=['x', 'y']) elif len(value.shape) == 1: - pass # Axes + pass # Axes else: failed = True else: From 5d5755b3cc3955d57256a5f016656d2733d21d46 Mon Sep 17 00:00:00 2001 From: Victor de Souza magalhaes Date: Wed, 26 Feb 2025 15:29:00 -0700 Subject: [PATCH 35/53] Added a chunk execution method --- src/astrohack/core/image_compare_tool.py | 32 ++++++++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/src/astrohack/core/image_compare_tool.py b/src/astrohack/core/image_compare_tool.py index 6859122d..13cf6393 100644 --- a/src/astrohack/core/image_compare_tool.py +++ b/src/astrohack/core/image_compare_tool.py @@ -296,3 +296,35 @@ def scatter_plot(self, destination, ref_image, dpi=300, display=False): ydata, f'{self.filename} [{self.unit}]', add_regression=True) close_figure(fig, 'Scatter plot against reference image', f'{destination}/{self.rootname}scatter.png', dpi, display) + + +def image_comparison_chunk(compare_params): + + image = FITSImage(compare_params['this_image'], compare_params['telescope_name']) + ref_image = FITSImage(compare_params['this_reference'], compare_params['telescope_name']) + plot_data = compare_params['plot_data'] + plot_percentuals = compare_params['plot_percentuals'] + plot_divided = compare_params['plot_divided_image'] + destination = compare_params['destination'] + colormap = compare_params['colormap'] + dpi = compare_params['dpi'] + display = compare_params['display'] + + if compare_params['comparison'] == 'direct': + image.compare_difference(ref_image) + image.plot_images(destination, plot_data, plot_percentuals, False, colormap=colormap, dpi=dpi, + display=display) + elif compare_params['comparison'] == 'scaled': + image.compare_scaled_difference(ref_image) + image.plot_images(destination, plot_data, plot_percentuals, plot_divided, colormap=colormap, dpi=dpi, + display=display) + else: + raise Exception(f'Unknown comparison type {compare_params["comparison"]}') + + if compare_params['export_to_fits']: + image.export_to_fits(destination) + + reference_node = xr.DataTree(name=ref_image.filename, data=ref_image.export_as_xds()) + tree_node = xr.DataTree(name=image.filename, data=image.export_as_xds(), children={'Reference': reference_node}) + + return tree_node From f7489b467254cccebf6836058b35a2acf6d662e7 Mon Sep 17 00:00:00 2001 From: Victor de Souza magalhaes Date: Wed, 26 Feb 2025 15:52:55 -0700 Subject: [PATCH 36/53] Added frontend version and improved chunk function --- src/astrohack/core/image_compare_tool.py | 14 +++++++++----- src/astrohack/image_comparison_tool.py | 24 ++++++++++++++++++++++++ 2 files changed, 33 insertions(+), 5 deletions(-) create mode 100644 src/astrohack/image_comparison_tool.py diff --git a/src/astrohack/core/image_compare_tool.py b/src/astrohack/core/image_compare_tool.py index 13cf6393..43808d41 100644 --- a/src/astrohack/core/image_compare_tool.py +++ b/src/astrohack/core/image_compare_tool.py @@ -63,7 +63,6 @@ def __init__(self, filename: str, telescope_name: str): def _init_as_fits(self, istokes=0, ichan=0): self.header, self.data = read_fits(self.filename, header_as_dict=True) - self.original_data = np.copy(self.data) self.fits_name = self.filename stokes_iaxis = get_stokes_axis_iaxis(self.header) @@ -74,12 +73,13 @@ def _init_as_fits(self, istokes=0, ichan=0): self.data = self.data[istokes, ichan, ...] else: self.data = self.data[ichan, istokes, ...] - elif len(self.data.shape) == 2: pass # image is already as expected else: raise Exception(f'FITS image has an unsupported shape: {self.data.shape}') + self.original_data = np.copy(self.data) + if 'AIPS' in self.header['ORIGIN']: self.x_axis, _, self.x_unit = get_axis_from_fits_header(self.header, 1, pixel_offset=False) self.y_axis, _, self.y_unit = get_axis_from_fits_header(self.header, 2, pixel_offset=False) @@ -301,7 +301,7 @@ def scatter_plot(self, destination, ref_image, dpi=300, display=False): def image_comparison_chunk(compare_params): image = FITSImage(compare_params['this_image'], compare_params['telescope_name']) - ref_image = FITSImage(compare_params['this_reference'], compare_params['telescope_name']) + ref_image = FITSImage(compare_params['this_reference_image'], compare_params['telescope_name']) plot_data = compare_params['plot_data'] plot_percentuals = compare_params['plot_percentuals'] plot_divided = compare_params['plot_divided_image'] @@ -324,7 +324,11 @@ def image_comparison_chunk(compare_params): if compare_params['export_to_fits']: image.export_to_fits(destination) - reference_node = xr.DataTree(name=ref_image.filename, data=ref_image.export_as_xds()) - tree_node = xr.DataTree(name=image.filename, data=image.export_as_xds(), children={'Reference': reference_node}) + if compare_params['plot_scatter']: + image.scatter_plot(destination, ref_image, dpi=dpi, display=display) + + img_node = xr.DataTree(name=image.filename, dataset=image.export_as_xds()) + ref_node = xr.DataTree(name=ref_image.filename, dataset=ref_image.export_as_xds()) + tree_node = xr.DataTree(name=image.rootname[:-1], children={'Reference': ref_node, 'Image': img_node}) return tree_node diff --git a/src/astrohack/image_comparison_tool.py b/src/astrohack/image_comparison_tool.py new file mode 100644 index 00000000..986804dd --- /dev/null +++ b/src/astrohack/image_comparison_tool.py @@ -0,0 +1,24 @@ +from astrohack.core.image_compare_tool import image_comparison_chunk + + +def compare_fits_image( + image: str, + reference_image: str, + telescope_name: str, + destination: str, + comparison: str = 'direct', + plot_data: bool = False, + plot_percentuals: bool = False, + plot_divided_image: bool = False, + plot_scatter: bool = True, + export_to_fits: bool = False, + colormap: str = 'viridis', + dpi: int = 300, + display: bool = False +): + param_dict = locals() + + param_dict['this_image'] = image + param_dict['this_reference_image'] = reference_image + + image_comparison_chunk(param_dict) From 3c484f1eed091e498d77a1e2d68eea28b24d66f2 Mon Sep 17 00:00:00 2001 From: Victor de Souza magalhaes Date: Wed, 26 Feb 2025 16:06:12 -0700 Subject: [PATCH 37/53] Some progress towards controlling the image comparison tool --- src/astrohack/image_comparison_tool.py | 28 ++++++++++++++++++++------ 1 file changed, 22 insertions(+), 6 deletions(-) diff --git a/src/astrohack/image_comparison_tool.py b/src/astrohack/image_comparison_tool.py index 986804dd..b424bf0c 100644 --- a/src/astrohack/image_comparison_tool.py +++ b/src/astrohack/image_comparison_tool.py @@ -1,9 +1,15 @@ +from typing import Union, List +import xarray as xr +import pathlib + +import toolviper.utils.logger as logger + from astrohack.core.image_compare_tool import image_comparison_chunk def compare_fits_image( - image: str, - reference_image: str, + image: Union[str, List[str]], + reference_image: Union[str, List[str]], telescope_name: str, destination: str, comparison: str = 'direct', @@ -16,9 +22,19 @@ def compare_fits_image( dpi: int = 300, display: bool = False ): - param_dict = locals() - param_dict['this_image'] = image - param_dict['this_reference_image'] = reference_image + if isinstance(image, str): + image = [image] + if isinstance(reference_image, str): + reference_image = [reference_image] + if len(image) != len(reference_image): + msg = 'List of reference images has a different size from the list of images' + logger.error(msg) + return - image_comparison_chunk(param_dict) + param_dict = locals() + pathlib.Path(param_dict['destination']).mkdir(exist_ok=True) + for i_img in range(len(image)): + param_dict['this_image'] = image[i_img] + param_dict['this_reference_image'] = reference_image[i_img] + image_comparison_chunk(param_dict) From 2b62f56c923bf97f86d9df5e2c634568c68c9805 Mon Sep 17 00:00:00 2001 From: Victor de Souza magalhaes Date: Wed, 26 Feb 2025 16:23:18 -0700 Subject: [PATCH 38/53] Added a new graph execution function --- src/astrohack/utils/graph.py | 33 +++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/src/astrohack/utils/graph.py b/src/astrohack/utils/graph.py index f989a19b..4f88f500 100644 --- a/src/astrohack/utils/graph.py +++ b/src/astrohack/utils/graph.py @@ -89,3 +89,36 @@ def compute_graph(looping_dict, chunk_function, param_dict, key_order, parallel= if parallel: dask.compute(delayed_list) return True + + +def compute_graph_from_lists(param_dict, chunk_function, looping_key_list, parallel=False): + """ + Creates and executes a graph based on entries in a parameter dictionary that are lists + Args: + param_dict: The parameter dictionary + chunk_function: The function for the operation chunk + looping_key_list: The keys that are lists in the parameter dictionaries over which to loop over + parallel: execute graph in parallel? + + Returns: + A list containing the returns of the calls to the chunk function. + """ + niter = len(param_dict[looping_key_list[0]]) + + delayed_list = [] + result_list = [] + for i_iter in range(niter): + this_param = param_dict.copy() + for key in looping_key_list: + this_param[f'this_{key}'] = param_dict[key][i_iter] + + if parallel: + delayed_list.append(dask.delayed(chunk_function)(dask.delayed(this_param))) + else: + delayed_list.append(0) + result_list.append(chunk_function(this_param)) + + if parallel: + result_list = dask.compute(delayed_list) + + return result_list From 70ac9c80e70ab7058e2f8f2264a931e465a8381d Mon Sep 17 00:00:00 2001 From: Victor de Souza magalhaes Date: Wed, 26 Feb 2025 16:42:25 -0700 Subject: [PATCH 39/53] Added creation of a datatree to disk. --- src/astrohack/image_comparison_tool.py | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/src/astrohack/image_comparison_tool.py b/src/astrohack/image_comparison_tool.py index b424bf0c..9ae4038b 100644 --- a/src/astrohack/image_comparison_tool.py +++ b/src/astrohack/image_comparison_tool.py @@ -5,6 +5,7 @@ import toolviper.utils.logger as logger from astrohack.core.image_compare_tool import image_comparison_chunk +from astrohack.utils.graph import compute_graph_from_lists def compare_fits_image( @@ -13,6 +14,7 @@ def compare_fits_image( telescope_name: str, destination: str, comparison: str = 'direct', + zarr_container_name: str = None, plot_data: bool = False, plot_percentuals: bool = False, plot_divided_image: bool = False, @@ -20,7 +22,8 @@ def compare_fits_image( export_to_fits: bool = False, colormap: str = 'viridis', dpi: int = 300, - display: bool = False + display: bool = False, + parallel: bool = False ): if isinstance(image, str): @@ -34,7 +37,13 @@ def compare_fits_image( param_dict = locals() pathlib.Path(param_dict['destination']).mkdir(exist_ok=True) - for i_img in range(len(image)): - param_dict['this_image'] = image[i_img] - param_dict['this_reference_image'] = reference_image[i_img] - image_comparison_chunk(param_dict) + + result_list = compute_graph_from_lists(param_dict, image_comparison_chunk, ['image', 'reference_image'], parallel) + + if zarr_container_name is not None: + root = xr.DataTree(name='Root') + for item in result_list: + tree_node = item[0] + root = root.assign({tree_node.name: tree_node}) + + root.to_zarr(zarr_container_name, mode='w', consolidated=True) From 5e4516caf370f17f95ac031e5448dfe677ab581e Mon Sep 17 00:00:00 2001 From: Victor de Souza magalhaes Date: Thu, 27 Feb 2025 09:59:18 -0700 Subject: [PATCH 40/53] Renaming --- .../core/{image_compare_tool.py => image_comparison_tool.py} | 0 src/astrohack/image_comparison_tool.py | 2 +- 2 files changed, 1 insertion(+), 1 deletion(-) rename src/astrohack/core/{image_compare_tool.py => image_comparison_tool.py} (100%) diff --git a/src/astrohack/core/image_compare_tool.py b/src/astrohack/core/image_comparison_tool.py similarity index 100% rename from src/astrohack/core/image_compare_tool.py rename to src/astrohack/core/image_comparison_tool.py diff --git a/src/astrohack/image_comparison_tool.py b/src/astrohack/image_comparison_tool.py index 9ae4038b..a91fc3ff 100644 --- a/src/astrohack/image_comparison_tool.py +++ b/src/astrohack/image_comparison_tool.py @@ -4,7 +4,7 @@ import toolviper.utils.logger as logger -from astrohack.core.image_compare_tool import image_comparison_chunk +from astrohack.core.image_comparison_tool import image_comparison_chunk from astrohack.utils.graph import compute_graph_from_lists From a6a7bab6c255f4fd9ce46af0d7403e6ce70906fa Mon Sep 17 00:00:00 2001 From: Victor de Souza magalhaes Date: Thu, 27 Feb 2025 10:18:22 -0700 Subject: [PATCH 41/53] Created class methods for proper initialization. --- src/astrohack/core/image_comparison_tool.py | 68 ++++++++++----------- 1 file changed, 34 insertions(+), 34 deletions(-) diff --git a/src/astrohack/core/image_comparison_tool.py b/src/astrohack/core/image_comparison_tool.py index 43808d41..1a88ee1b 100644 --- a/src/astrohack/core/image_comparison_tool.py +++ b/src/astrohack/core/image_comparison_tool.py @@ -22,48 +22,52 @@ def test_image(fits_image): class FITSImage: - def __init__(self, filename: str, telescope_name: str): + def __init__(self): """ - Initializes the FITSImage object from a file on disk - Args: - filename: Name of the file on disk, may be .FITS file or a .zarr xds with a disk representation of a \ - FITSImage object - telescope_name: Name of the telescope used on the images so that masking can be properly applied. + Blank slate initialization of the FITSImage object """ - # Initialization from parameters - self.filename = filename - self.telescope_name = telescope_name - self.rootname = '.'.join(filename.split('.')[:-1])+'.' - - # Blank slate initialization - self.header = None - self.data = None + # Attributes: + self.filename = None + self.telescope_name = None + self.rootname = None self.factor = 1.0 - self.residuals = None - self.residuals_percent = None - self.divided_image = None self.reference_name = None self.resampled = False + + # Metadata + self.header = None + self.unit = None self.x_axis = None self.y_axis = None self.original_x_axis = None self.original_y_axis = None self.x_unit = None self.y_unit = None - self.unit = None - self.fits_name = None + + # Data variables self.original_data = None + self.data = None + self.residuals = None + self.residuals_percent = None + self.divided_image = None - if '.FITS' in filename.upper(): - self._init_as_fits() - elif '.zarr' in filename: - self._init_as_xds() - else: - raise Exception(f"Don't know how to read {filename}") + @classmethod + def from_xds(cls, xds): + return_obj = cls() + return_obj._init_as_xds(xds) + return return_obj + + @classmethod + def from_fits_file(cls, fits_filename, telescope_name): + return_obj = cls() + return_obj._init_as_fits(fits_filename, telescope_name) + return return_obj - def _init_as_fits(self, istokes=0, ichan=0): + def _init_as_fits(self, fits_filename, telescope_name, istokes=0, ichan=0): + self.filename = fits_filename + self.telescope_name = telescope_name + self.rootname = '.'.join(fits_filename.split('.')[:-1])+'.' self.header, self.data = read_fits(self.filename, header_as_dict=True) - self.fits_name = self.filename stokes_iaxis = get_stokes_axis_iaxis(self.header) self.unit = self.header['BUNIT'] @@ -95,9 +99,7 @@ def _init_as_fits(self, istokes=0, ichan=0): self.original_x_axis = np.copy(self.x_axis) self.original_y_axis = np.copy(self.y_axis) - def _init_as_xds(self): - filename = self.filename - xds = xr.open_zarr(self.filename) + def _init_as_xds(self, xds): for key in xds.attrs: setattr(self, key, xds.attrs[key]) @@ -109,8 +111,6 @@ def _init_as_xds(self): for key, value in xds.items(): setattr(self, str(key), xds[key].values) - self.filename = filename - def _create_base_mask(self): telescope_obj = Telescope(self.telescope_name) self.base_mask = create_aperture_mask(self.x_axis, self.y_axis, telescope_obj.inlim, telescope_obj.oulim, @@ -300,8 +300,8 @@ def scatter_plot(self, destination, ref_image, dpi=300, display=False): def image_comparison_chunk(compare_params): - image = FITSImage(compare_params['this_image'], compare_params['telescope_name']) - ref_image = FITSImage(compare_params['this_reference_image'], compare_params['telescope_name']) + image = FITSImage.from_fits_file(compare_params['this_image'], compare_params['telescope_name']) + ref_image = FITSImage.from_fits_file(compare_params['this_reference_image'], compare_params['telescope_name']) plot_data = compare_params['plot_data'] plot_percentuals = compare_params['plot_percentuals'] plot_divided = compare_params['plot_divided_image'] From 76e88d64bac57dc08953baff4a273cc769736d75 Mon Sep 17 00:00:00 2001 From: Victor de Souza magalhaes Date: Thu, 27 Feb 2025 11:03:08 -0700 Subject: [PATCH 42/53] Documentation for FITSImage class --- src/astrohack/core/image_comparison_tool.py | 164 ++++++++++++++++++++ src/astrohack/image_comparison_tool.py | 2 +- 2 files changed, 165 insertions(+), 1 deletion(-) diff --git a/src/astrohack/core/image_comparison_tool.py b/src/astrohack/core/image_comparison_tool.py index 1a88ee1b..e9b8b0c0 100644 --- a/src/astrohack/core/image_comparison_tool.py +++ b/src/astrohack/core/image_comparison_tool.py @@ -53,17 +53,60 @@ def __init__(self): @classmethod def from_xds(cls, xds): + """ + Initialize a FITSImage object using as a base a Xarray dataset + Args: + xds: Xarray dataset + + Returns: + FITSImage object initialized from a xds + """ return_obj = cls() return_obj._init_as_xds(xds) return return_obj @classmethod def from_fits_file(cls, fits_filename, telescope_name): + """ + Initialize a FITSImage object using as a base a FITS file. + Args: + fits_filename: FITS file on disk + telescope_name: Name of the telescope used + + Returns: + FITSImage object initialized from a FITS file + """ return_obj = cls() return_obj._init_as_fits(fits_filename, telescope_name) return return_obj + @classmethod + def from_zarr(cls, zarr_filename): + """ + Initialize a FITSImage object using as a base a Xarray dataset store on disk in a zarr container + Args: + zarr_filename: Xarray dataset on disk as a zarr container + + Returns: + FITSImage object initialized from a xds + """ + return_obj = cls() + xds = xr.open_zarr(zarr_filename) + return_obj._init_as_xds(xds) + return return_obj + def _init_as_fits(self, fits_filename, telescope_name, istokes=0, ichan=0): + """ + Backend for FITSImage.from_fits_file + Args: + fits_filename: FITS file on disk + telescope_name: Name of the telescope used + istokes: Stokes axis element to be fetched, should always be zero (singleton stokes axis or fetching I) + ichan: Channel axis element to be fetched, should be zero for most cases, unless image has multiple channels + + Returns: + None + """ self.filename = fits_filename self.telescope_name = telescope_name self.rootname = '.'.join(fits_filename.split('.')[:-1])+'.' @@ -100,6 +143,13 @@ def _init_as_fits(self, fits_filename, telescope_name, istokes=0, ichan=0): self.original_y_axis = np.copy(self.y_axis) def _init_as_xds(self, xds): + """ + Backend for FITSImage.from_xds + Args: + xds: Xarray DataSet + Returns: + None + """ for key in xds.attrs: setattr(self, key, xds.attrs[key]) @@ -112,12 +162,25 @@ def _init_as_xds(self, xds): setattr(self, str(key), xds[key].values) def _create_base_mask(self): + """ + Create a base mask based on telescope parameters such as arm shadows. + Returns: + None + """ telescope_obj = Telescope(self.telescope_name) self.base_mask = create_aperture_mask(self.x_axis, self.y_axis, telescope_obj.inlim, telescope_obj.oulim, arm_width=telescope_obj.arm_shadow_width, arm_angle=telescope_obj.arm_shadow_rotation) def resample(self, ref_image): + """ + Resamples the data on this object onto the grid in ref_image + Args: + ref_image: Reference FITSImage object + + Returns: + None + """ test_image(ref_image) x_mesh_orig, y_mesh_orig = np.meshgrid(self.x_axis, self.y_axis, indexing='ij') x_mesh_dest, y_mesh_dest = np.meshgrid(ref_image.x_axis, ref_image.y_axis, indexing='ij') @@ -132,6 +195,14 @@ def resample(self, ref_image): self.resampled = True def compare_difference(self, ref_image): + """ + Does the difference comparison between self and ref_image. + Args: + ref_image: Reference FITSImage object + + Returns: + None + """ test_image(ref_image) if not self.image_has_same_sampling(ref_image): self.resample(ref_image) @@ -141,6 +212,15 @@ def compare_difference(self, ref_image): self.reference_name = ref_image.filename def compare_scaled_difference(self, ref_image, rejection=10): + """ + Does the scaled difference comparison between self and ref_image. + Args: + ref_image: Reference FITSImage object + rejection: rejection level for scaling factor + + Returns: + None + """ test_image(ref_image) if not self.image_has_same_sampling(ref_image): self.resample(ref_image) @@ -151,14 +231,44 @@ def compare_scaled_difference(self, ref_image, rejection=10): self.compare_difference(ref_image) def image_has_same_sampling(self, ref_image): + """ + Tests if self has the same X and Y sampling as ref_image + Args: + ref_image: Reference FITSImage object + + Returns: + True or False + """ test_image(ref_image) return are_axes_equal(self.x_axis, ref_image.x_axis) and are_axes_equal(self.y_axis, ref_image.y_axis) def _mask_array(self, image_array): + """ + Applies base mask to image_array + Args: + image_array: Data array to be masked + + Returns: + Masked array + """ return np.where(self.base_mask, image_array, np.nan) def plot_images(self, destination, plot_data=False, plot_percentuals=False, plot_divided_image=False, colormap='viridis', dpi=300, display=False): + """ + Plot image contents of the FITSImage object, always plots the residuals when called + Args: + destination: Location onto which save plot files + plot_data: Also plot data array? + plot_percentuals: Also plot percentual residuals array? + plot_divided_image: Also plot divided image? + colormap: Colormap name for image plots + dpi: png resolution on disk + display: Show interactive view of plots + + Returns: + None + """ extent = compute_extent(self.x_axis, self.y_axis, 0.0) cmap = get_proper_color_map(colormap) @@ -190,6 +300,23 @@ def plot_images(self, destination, plot_data=False, plot_percentuals=False, dpi, display, add_statistics=True) def _plot_map(self, data, title, zlabel, filename, cmap, extent, zscale, dpi, display, add_statistics=False): + """ + Backend for plot_images + Args: + data: Data array to be plotted + title: Title to appear on plot + zlabel: Label for the colorbar + filename: name for the png file on disk + cmap: Colormap object for plots + extent: extents of the X and Y axes + zscale: Constraints on the Z axes. + dpi: png resolution on disk + display: Show interactive view of plots + add_statistics: Add simple statistics to plot's subtitle + + Returns: + None + """ fig, ax = plt.subplots(1, 1, figsize=[10, 8]) if zscale == 'symmetrical': scale = max(np.abs(np.nanmin(data)), np.abs(np.nanmax(data))) @@ -212,6 +339,11 @@ def _plot_map(self, data, title, zlabel, filename, cmap, extent, zscale, dpi, di close_figure(fig, title, filename, dpi, display) def export_as_xds(self): + """ + Create a Xarray DataSet from the FITSImage object + Returns: + Xarray DataSet + """ xds = xr.Dataset() obj_dict = vars(self) @@ -239,10 +371,23 @@ def export_as_xds(self): return xds def to_zarr(self, zarr_filename): + """ + Saves a xds representation of self on disk using the zarr format. + Args: + zarr_filename: Name for the zarr container on disk + + Returns: + None + """ xds = self.export_as_xds() xds.to_zarr(zarr_filename, mode="w", compute=True, consolidated=True) def __repr__(self): + """ + Print method + Returns: + A String summary of the current status of self. + """ obj_dict = vars(self) outstr = '' for key, value in obj_dict.items(): @@ -256,6 +401,14 @@ def __repr__(self): return outstr def export_to_fits(self, destination): + """ + Export internal images to FITS files. + Args: + destination: location to store FITS files + + Returns: + None + """ pathlib.Path(destination).mkdir(exist_ok=True) ext_fits = '.fits' out_header = self.header.copy() @@ -281,6 +434,17 @@ def export_to_fits(self, destination): write_fits(out_header, key, np.fliplr(value.astype(float)), filename, unit, reorder_axis=False) def scatter_plot(self, destination, ref_image, dpi=300, display=False): + """ + Produce a scatter plot of self.data agains ref_image.data + Args: + destination: Location to store scatter plot + ref_image: Reference FITSImage object + dpi: png resolution on disk + display: Show interactive view of plot + + Returns: + None + """ test_image(ref_image) if not self.image_has_same_sampling(ref_image): self.resample(ref_image) diff --git a/src/astrohack/image_comparison_tool.py b/src/astrohack/image_comparison_tool.py index a91fc3ff..d4673dd4 100644 --- a/src/astrohack/image_comparison_tool.py +++ b/src/astrohack/image_comparison_tool.py @@ -8,7 +8,7 @@ from astrohack.utils.graph import compute_graph_from_lists -def compare_fits_image( +def compare_fits_images( image: Union[str, List[str]], reference_image: Union[str, List[str]], telescope_name: str, From b04278afe6f074cdc0643cb982ac2a40b832e540 Mon Sep 17 00:00:00 2001 From: Victor de Souza magalhaes Date: Thu, 27 Feb 2025 11:47:51 -0700 Subject: [PATCH 43/53] Further backend documentation. --- src/astrohack/core/image_comparison_tool.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/src/astrohack/core/image_comparison_tool.py b/src/astrohack/core/image_comparison_tool.py index e9b8b0c0..0fbe73c1 100644 --- a/src/astrohack/core/image_comparison_tool.py +++ b/src/astrohack/core/image_comparison_tool.py @@ -295,7 +295,8 @@ def plot_images(self, destination, plot_data=False, plot_percentuals=False, if plot_divided_image: if self.divided_image is None: raise Exception("Cannot plot a divided image that does not exist.") - self._plot_map(self._mask_array(self.divided_image), f'Divided image, ref={self.reference_name}', + self._plot_map(self._mask_array(self.divided_image), + f'Divided image, ref={self.reference_name}, scaling={self.factor:.4f}', f'Division [ ]', f'{base_name}divided.png', cmap, extent, [None, None], dpi, display, add_statistics=True) @@ -463,6 +464,14 @@ def scatter_plot(self, destination, ref_image, dpi=300, display=False): def image_comparison_chunk(compare_params): + """ + Chunk function for parallel execution of the image comparison tool. + Args: + compare_params: Parameter dictionary for workflow control. + + Returns: + A DataTree containing the Image and its reference Image. + """ image = FITSImage.from_fits_file(compare_params['this_image'], compare_params['telescope_name']) ref_image = FITSImage.from_fits_file(compare_params['this_reference_image'], compare_params['telescope_name']) From a6ffb15dab061fa226bce5237ebef93d502b26ba Mon Sep 17 00:00:00 2001 From: Victor de Souza magalhaes Date: Thu, 27 Feb 2025 11:48:32 -0700 Subject: [PATCH 44/53] Frontend documentation and added DataTree return. --- src/astrohack/image_comparison_tool.py | 97 ++++++++++++++++++++++++-- 1 file changed, 92 insertions(+), 5 deletions(-) diff --git a/src/astrohack/image_comparison_tool.py b/src/astrohack/image_comparison_tool.py index d4673dd4..90b853c1 100644 --- a/src/astrohack/image_comparison_tool.py +++ b/src/astrohack/image_comparison_tool.py @@ -25,6 +25,91 @@ def compare_fits_images( display: bool = False, parallel: bool = False ): + """ + Compares a set of images to a set of reference images. + + :param image: FITS image or list of FITS images to be compared. + :type image: list or str + + :param reference_image: FITS image or list of FITS images that serve as references. + :type reference_image: list or str + + :param telescope_name: Name of the telescope used. Used for masking. + :type telescope_name: str + + :param destination: Name of directory onto which save plots + :type destination: str + + :param comparison: Type of comparison to be made between images, "direct" or "scaled", default is "direct". + :type comparison: str, optional + + :param zarr_container_name: Name of the Zarr container to contain the created datatree, default is None, i.e. \ + DataTree is not saved to disk. + :type zarr_container_name: str, optional + + :param plot_data: Plot the data array used in the comparison, default is False. + :type plot_data: bool, optional + + :param plot_percentuals: Plot the residuals in percent of reference image as well, default is False. + :type plot_percentuals: bool, optional + + :param plot_divided_image: Plot the divided image between Image and its reference, default is False. + :type plot_divided_image: bool, optional + + :param plot_scatter: Make a scatter plot of the Image against its reference image, default is True. + :type plot_scatter: bool, optional + + :param export_to_fits: Export created images to FITS files inside destination, default is False. + :type export_to_fits: bool, optional + + :param colormap: Colormap to be used on image plots, default is "viridis". + :type colormap: str, optional + + :param dpi: dots per inch to be used in plots, default is 300. + :type dpi: int, optional + + :param display: Display plots inline or suppress, defaults to True + :type display: bool, optional + + :param parallel: If True will use an existing astrohack client to do comparison in parallel, default is False + :type parallel: bool, optional + + :return: DataTree object containing all the comparisons executed + :rtype: xr.DataTree + + .. _Description: + Compares pairs of FITS images pixel by pixel using a mask based on telescope parameters to exclude problematic \ + regions such as shadows caused by the secondary mirror or the arms supporting it. By default, 2 products are \ + produced, a plot of the residuals image, i.e. (Reference - Image) and a scatter plot of the Reference against the \ + Image. If necessary a resample of Image is conducted to allow for pixel by pixel comparison. + + .. rubric:: Comparison: + Two types of comparison between the images are available: + - *direct*: Where the residuals are simply computed as Reference - Image. + - *scaled*: Where the residuals are Reference - Factor * Image, with Factor = median(Reference/Image). + + .. rubric:: Plots: + A plot of the residuals of the comparison is always produced. + However, a few extra plots can be produced and their production is controlled by the plot_ parameters, these are: + - *plot_data*: Activates plotting of the data used in the comparison, default is False as this is the data on the \ + FITS file. + - *plot_percentuals*: Activates the plotting of the residuals as a perdentage of the Reference Image, default is \ + False as this is just another view on the residuals. + -*plot_divided_image*: Activates the plotting of Reference/Image, default is False. This plot is only available \ + when using "scaled" comparison. + -*plot_scatter*: Activates the creation of a scatter plot of Reference vs Image, with a linear regression, default \ + is True. + + .. rubric:: Storage on disk: + By default, this function only produces plots, but this can be changed using two parameters: + - *zarr_container_name*: If this parameter is not None a Zarr container will be created on disk with the contents \ + of the produced DataTree. + - *export_to_fits*: If set to True will produce FITS files of the produced images and store them at *destination*. + + .. rubric:: Return type: + This funtion returns a Xarray DataTree containing the Xarray DataSets that represent Image and Reference. The nodes \ + in this DataTree are labelled according to the filenames given as input for easier navigation. + """ if isinstance(image, str): image = [image] @@ -40,10 +125,12 @@ def compare_fits_images( result_list = compute_graph_from_lists(param_dict, image_comparison_chunk, ['image', 'reference_image'], parallel) - if zarr_container_name is not None: - root = xr.DataTree(name='Root') - for item in result_list: - tree_node = item[0] - root = root.assign({tree_node.name: tree_node}) + root = xr.DataTree(name='Root') + for item in result_list: + tree_node = item[0] + root = root.assign({tree_node.name: tree_node}) + if zarr_container_name is not None: root.to_zarr(zarr_container_name, mode='w', consolidated=True) + + return root From c94d4e78533f338b3bb348c12b49a1739d690a84 Mon Sep 17 00:00:00 2001 From: Victor de Souza magalhaes Date: Thu, 27 Feb 2025 11:50:44 -0700 Subject: [PATCH 45/53] Added image comparison tool to defautl astrohack import. --- src/astrohack/__init__.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/astrohack/__init__.py b/src/astrohack/__init__.py index ed095d08..979771c6 100644 --- a/src/astrohack/__init__.py +++ b/src/astrohack/__init__.py @@ -14,6 +14,7 @@ from .locit import * from .extract_locit import * from .cassegrain_ray_tracing import * +from .image_comparison_tool import * # This installs a slick, informational tracebacks logger from rich.traceback import install From 5786bcf06981416d6c31b617e204d627cbad176e Mon Sep 17 00:00:00 2001 From: Victor de Souza magalhaes Date: Thu, 27 Feb 2025 11:50:57 -0700 Subject: [PATCH 46/53] Added image comparison tool to API. --- docs/api.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/api.rst b/docs/api.rst index 42c391e6..2b2800d1 100755 --- a/docs/api.rst +++ b/docs/api.rst @@ -16,3 +16,4 @@ Under development. _api/autoapi/astrohack/extract_locit/index _api/autoapi/astrohack/locit/index _api/autoapi/astrohack/cassegrain_ray_tracing/index + _api/autoapi/astrohack/image_comparison_tool/index From 742ecf53d626986ac18183da8c8b00ac5b124238 Mon Sep 17 00:00:00 2001 From: Victor de Souza magalhaes Date: Thu, 27 Feb 2025 11:56:01 -0700 Subject: [PATCH 47/53] Changes in format for documentation --- src/astrohack/image_comparison_tool.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/astrohack/image_comparison_tool.py b/src/astrohack/image_comparison_tool.py index 90b853c1..214b8439 100644 --- a/src/astrohack/image_comparison_tool.py +++ b/src/astrohack/image_comparison_tool.py @@ -85,26 +85,26 @@ def compare_fits_images( .. rubric:: Comparison: Two types of comparison between the images are available: - - *direct*: Where the residuals are simply computed as Reference - Image. - - *scaled*: Where the residuals are Reference - Factor * Image, with Factor = median(Reference/Image). + * *direct*: Where the residuals are simply computed as Reference - Image. + * *scaled*: Where the residuals are Reference - Factor * Image, with Factor = median(Reference/Image). .. rubric:: Plots: A plot of the residuals of the comparison is always produced. However, a few extra plots can be produced and their production is controlled by the plot_ parameters, these are: - - *plot_data*: Activates plotting of the data used in the comparison, default is False as this is the data on the \ + * *plot_data*: Activates plotting of the data used in the comparison, default is False as this is the data on the \ FITS file. - - *plot_percentuals*: Activates the plotting of the residuals as a perdentage of the Reference Image, default is \ + * *plot_percentuals*: Activates the plotting of the residuals as a perdentage of the Reference Image, default is \ False as this is just another view on the residuals. - -*plot_divided_image*: Activates the plotting of Reference/Image, default is False. This plot is only available \ + * *plot_divided_image*: Activates the plotting of Reference/Image, default is False. This plot is only available \ when using "scaled" comparison. - -*plot_scatter*: Activates the creation of a scatter plot of Reference vs Image, with a linear regression, default \ + * *plot_scatter*: Activates the creation of a scatter plot of Reference vs Image, with a linear regression, default \ is True. .. rubric:: Storage on disk: By default, this function only produces plots, but this can be changed using two parameters: - - *zarr_container_name*: If this parameter is not None a Zarr container will be created on disk with the contents \ + * *zarr_container_name*: If this parameter is not None a Zarr container will be created on disk with the contents \ of the produced DataTree. - - *export_to_fits*: If set to True will produce FITS files of the produced images and store them at *destination*. + * *export_to_fits*: If set to True will produce FITS files of the produced images and store them at *destination*. .. rubric:: Return type: This funtion returns a Xarray DataTree containing the Xarray DataSets that represent Image and Reference. The nodes \ From 4ebbf68213fdb61c595e4586e8517e9950345157 Mon Sep 17 00:00:00 2001 From: Victor de Souza magalhaes Date: Thu, 27 Feb 2025 12:01:29 -0700 Subject: [PATCH 48/53] Changes in format for documentation --- src/astrohack/image_comparison_tool.py | 1 + src/astrohack/panel.py | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/src/astrohack/image_comparison_tool.py b/src/astrohack/image_comparison_tool.py index 214b8439..5aebcbc0 100644 --- a/src/astrohack/image_comparison_tool.py +++ b/src/astrohack/image_comparison_tool.py @@ -85,6 +85,7 @@ def compare_fits_images( .. rubric:: Comparison: Two types of comparison between the images are available: + * *direct*: Where the residuals are simply computed as Reference - Image. * *scaled*: Where the residuals are Reference - Factor * Image, with Factor = median(Reference/Image). diff --git a/src/astrohack/panel.py b/src/astrohack/panel.py index 88a3e1c1..8c8d1445 100644 --- a/src/astrohack/panel.py +++ b/src/astrohack/panel.py @@ -52,7 +52,7 @@ def panel( passing a dictionary, default is 3 (appropriate for sigma clipping) :type clip_level: float, dict, optional - :param exclude_shadows: Exclude regions with significant shadowing from analysis, e.g. secondary supporting arms, + :param exclude_shadows: Exclude regions with significant shadowing from analysis, e.g. secondary supporting arms, \ default is True. :type exclude_shadows: bool, optional From 4f1c69285e85c14c63197fee2b050aa5eef000f8 Mon Sep 17 00:00:00 2001 From: Victor de Souza magalhaes Date: Thu, 27 Feb 2025 12:20:23 -0700 Subject: [PATCH 49/53] Changes in format for documentation --- src/astrohack/image_comparison_tool.py | 28 +++++++++++++------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/src/astrohack/image_comparison_tool.py b/src/astrohack/image_comparison_tool.py index 5aebcbc0..9674e7d4 100644 --- a/src/astrohack/image_comparison_tool.py +++ b/src/astrohack/image_comparison_tool.py @@ -85,27 +85,27 @@ def compare_fits_images( .. rubric:: Comparison: Two types of comparison between the images are available: - - * *direct*: Where the residuals are simply computed as Reference - Image. - * *scaled*: Where the residuals are Reference - Factor * Image, with Factor = median(Reference/Image). + - *direct*: Where the residuals are simply computed as Reference - Image. + - *scaled*: Where the residuals are Reference - Factor * Image, with Factor = median(Reference/Image). .. rubric:: Plots: A plot of the residuals of the comparison is always produced. However, a few extra plots can be produced and their production is controlled by the plot_ parameters, these are: - * *plot_data*: Activates plotting of the data used in the comparison, default is False as this is the data on the \ - FITS file. - * *plot_percentuals*: Activates the plotting of the residuals as a perdentage of the Reference Image, default is \ - False as this is just another view on the residuals. - * *plot_divided_image*: Activates the plotting of Reference/Image, default is False. This plot is only available \ - when using "scaled" comparison. - * *plot_scatter*: Activates the creation of a scatter plot of Reference vs Image, with a linear regression, default \ - is True. + - *plot_data*: Activates plotting of the data used in the comparison, default is False as this is the data on \ + the FITS file. + - *plot_percentuals*: Activates the plotting of the residuals as a perdentage of the Reference Image, default \ + is False as this is just another view on the residuals. + - *plot_divided_image*: Activates the plotting of Reference/Image, default is False. This plot is only \ + available when using "scaled" comparison. + - *plot_scatter*: Activates the creation of a scatter plot of Reference vs Image, with a linear regression, \ + default is True. .. rubric:: Storage on disk: By default, this function only produces plots, but this can be changed using two parameters: - * *zarr_container_name*: If this parameter is not None a Zarr container will be created on disk with the contents \ - of the produced DataTree. - * *export_to_fits*: If set to True will produce FITS files of the produced images and store them at *destination*. + - *zarr_container_name*: If this parameter is not None a Zarr container will be created on disk with the \ + contents of the produced DataTree. + - *export_to_fits*: If set to True will produce FITS files of the produced images and store them at \ + *destination*. .. rubric:: Return type: This funtion returns a Xarray DataTree containing the Xarray DataSets that represent Image and Reference. The nodes \ From 8d04fee042aa6f3e678080d0765cf72e71424ea2 Mon Sep 17 00:00:00 2001 From: Victor de Souza magalhaes Date: Thu, 27 Feb 2025 12:27:26 -0700 Subject: [PATCH 50/53] Changes in format for documentation --- src/astrohack/image_comparison_tool.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/astrohack/image_comparison_tool.py b/src/astrohack/image_comparison_tool.py index 9674e7d4..11844832 100644 --- a/src/astrohack/image_comparison_tool.py +++ b/src/astrohack/image_comparison_tool.py @@ -90,7 +90,7 @@ def compare_fits_images( .. rubric:: Plots: A plot of the residuals of the comparison is always produced. - However, a few extra plots can be produced and their production is controlled by the plot_ parameters, these are: + However, a few extra plots can be produced and their production is controlled by the *plot_* parameters, these are: - *plot_data*: Activates plotting of the data used in the comparison, default is False as this is the data on \ the FITS file. - *plot_percentuals*: Activates the plotting of the residuals as a perdentage of the Reference Image, default \ From 23a9403e2faf05ccecc9574d8603139a993c274b Mon Sep 17 00:00:00 2001 From: Victor de Souza magalhaes Date: Thu, 27 Feb 2025 13:53:52 -0700 Subject: [PATCH 51/53] Added parameter checking for image comparison tool --- .../config/image_comparison_tool.param.json | 85 +++++++++++++++++++ src/astrohack/image_comparison_tool.py | 5 ++ 2 files changed, 90 insertions(+) create mode 100644 src/astrohack/config/image_comparison_tool.param.json diff --git a/src/astrohack/config/image_comparison_tool.param.json b/src/astrohack/config/image_comparison_tool.param.json new file mode 100644 index 00000000..49b7e5b4 --- /dev/null +++ b/src/astrohack/config/image_comparison_tool.param.json @@ -0,0 +1,85 @@ +{ + "compare_fits_images":{ + "image":{ + "nullable": false, + "required": true, + "struct_type": ["str"], + "type": ["string", "list"] + }, + "reference_image":{ + "nullable": false, + "required": true, + "struct_type": ["str"], + "type": ["string", "list"] + }, + "telescope_name":{ + "nullable": false, + "required": true, + "type": ["string"] + }, + "destination":{ + "nullable": false, + "required": true, + "type": ["string"] + }, + "comparison":{ + "allowed": ["direct", "scaled"], + "nullable": false, + "required": false, + "type": ["string"] + }, + "zarr_container_name":{ + "nullable": true, + "required": false, + "type": ["string"] + }, + "plot_data":{ + "nullable": false, + "required": false, + "type": ["boolean"] + }, + "plot_percentuals":{ + "nullable": false, + "required": false, + "type": ["boolean"] + }, + "plot_divided_image":{ + "nullable": false, + "required": false, + "type": ["boolean"] + }, + "plot_scatter":{ + "nullable": false, + "required": false, + "type": ["boolean"] + }, + "export_to_fits":{ + "nullable": false, + "required": false, + "type": ["boolean"] + }, + "colormap":{ + "nullable": false, + "required": false, + "type": ["string"], + "check allowed with": "colormaps" + }, + "dpi":{ + "nullable": false, + "required": false, + "type": ["int"], + "min": 1, + "max": 1200 + }, + "display":{ + "nullable": false, + "required": false, + "type": ["boolean"] + }, + "parallel":{ + "nullable": false, + "required": false, + "type": ["boolean"] + } + } +} diff --git a/src/astrohack/image_comparison_tool.py b/src/astrohack/image_comparison_tool.py index 11844832..ab8e3351 100644 --- a/src/astrohack/image_comparison_tool.py +++ b/src/astrohack/image_comparison_tool.py @@ -3,11 +3,16 @@ import pathlib import toolviper.utils.logger as logger +import toolviper from astrohack.core.image_comparison_tool import image_comparison_chunk from astrohack.utils.graph import compute_graph_from_lists +from astrohack.utils.validation import custom_plots_checker +@toolviper.utils.parameter.validate( + custom_checker=custom_plots_checker +) def compare_fits_images( image: Union[str, List[str]], reference_image: Union[str, List[str]], From d4820e6e5dc8b76647f278aebcc207222b79b9e4 Mon Sep 17 00:00:00 2001 From: Victor de Souza magalhaes Date: Thu, 27 Feb 2025 14:03:36 -0700 Subject: [PATCH 52/53] Updated reference values in test_panel.py and test_class_antenna_surface.py --- tests/unit/test_class_antenna_surface.py | 9 +++++---- tests/unit/test_panel.py | 2 +- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/tests/unit/test_class_antenna_surface.py b/tests/unit/test_class_antenna_surface.py index 2465ef16..a59ddcb6 100644 --- a/tests/unit/test_class_antenna_surface.py +++ b/tests/unit/test_class_antenna_surface.py @@ -90,8 +90,8 @@ def test_compile_panel_points_ringed(self): """ Tests that a point falls into the correct panel and that this panel has the correct number of samples """ - compvaluep0 = [3.2456341911764706, 0.7755055147058822, 198, 268, 0.00045656206805518506] - compnsampp0 = 120 + compvaluep0 = [3.3030790441176467, 0.43083639705882354, 197, 262, 0.00025600619999005243] + compnsampp0 = 179 self.tant.compile_panel_points() assert len(self.tant.panels[0].samples) == compnsampp0, 'Number of samples in panel is different from reference' assert self.tant.panels[0].samples[0] == PanelPoint(*compvaluep0), ('Point data in Panel is different from what' @@ -101,9 +101,10 @@ def test_fit_surface(self): """ Tests that fitting results for two panels match the reference """ - solveparsp0 = [0.00035746, 0.00020089, -0.0008455 ] - solveparsp30 = [ 0.00039911, -0.00041468, -0.0007079] + solveparsp0 = [0.00032385, 0.00037302, -0.00092492] + solveparsp30 = [0.00038098, -0.00039892, -0.00067244] self.tant.fit_surface() + assert len(self.tant.panels[0].model.parameters) == len(solveparsp0), ('Fitted results have a different length' ' from reference') for i in range(len(solveparsp30)): diff --git a/tests/unit/test_panel.py b/tests/unit/test_panel.py index 0f63952b..b4ea8f9f 100644 --- a/tests/unit/test_panel.py +++ b/tests/unit/test_panel.py @@ -202,7 +202,7 @@ def test_panel_absolute_clip(self): telescope = Telescope('vla') radius = panel_mds["ant_ea25"]["ddi_0"]['RADIUS'].values - dish_mask = np.where(radius < telescope.diam/2, 1.0, 0) + dish_mask = np.where(radius < telescope.oulim, 1.0, 0) dish_mask = np.where(radius < telescope.inlim, 0, dish_mask) nvalid_pix = np.sum(dish_mask) From abe134dc95b24051101e8acae76a45528b4f9440 Mon Sep 17 00:00:00 2001 From: Victor de Souza magalhaes Date: Thu, 27 Feb 2025 14:25:26 -0700 Subject: [PATCH 53/53] Updated reference values in test_stakeholder_vla.py --- tests/stakeholder/test_stakeholder_vla.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/tests/stakeholder/test_stakeholder_vla.py b/tests/stakeholder/test_stakeholder_vla.py index ae689d22..82346858 100644 --- a/tests/stakeholder/test_stakeholder_vla.py +++ b/tests/stakeholder/test_stakeholder_vla.py @@ -242,5 +242,9 @@ def test_holography_pipeline(set_data): exclude_shadows=False ) - assert verify_panel_shifts(data_dir=str(set_data), ref_mean_shift=reference_dict["vla"]["offsets"]), "Verify panel shifts" - #assert verify_panel_shifts(data_dir=str(set_data)), "Verify panel shifts" + reference_shifts = np.array([-91.6455227, 61.69666059, 4.39843319, 122.26547831]) + assert verify_panel_shifts(data_dir=str(set_data), ref_mean_shift=reference_shifts), "Verify panel shifts" + # This test using reference values is very hard to be updated, using this hardcoded reference_shifts is a + # temporary work around + # assert verify_panel_shifts(data_dir=str(set_data), ref_mean_shift=reference_dict["vla"]["offsets"]), \ + # "Verify panel shifts"