From f115a522a5af543da296a021d34b749ab6eacb53 Mon Sep 17 00:00:00 2001 From: Jeff Jennings Date: Wed, 29 Nov 2023 17:31:57 -0500 Subject: [PATCH 1/7] new script: input_output --- src/mpol/input_output.py | 4 ++++ 1 file changed, 4 insertions(+) create mode 100644 src/mpol/input_output.py diff --git a/src/mpol/input_output.py b/src/mpol/input_output.py new file mode 100644 index 00000000..36f08f9d --- /dev/null +++ b/src/mpol/input_output.py @@ -0,0 +1,4 @@ +import numpy as np + +from astropy.io import fits + From c0dcc1bc35045209f896740f5f0aeacb48801311 Mon Sep 17 00:00:00 2001 From: Jeff Jennings Date: Wed, 29 Nov 2023 17:32:11 -0500 Subject: [PATCH 2/7] input_output: new class: ProcessFitsImage --- src/mpol/input_output.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/src/mpol/input_output.py b/src/mpol/input_output.py index 36f08f9d..45093d4e 100644 --- a/src/mpol/input_output.py +++ b/src/mpol/input_output.py @@ -2,3 +2,19 @@ from astropy.io import fits +class ProcessFitsImage: + """ + Utilities for loading and retrieving metrics of a .fits image + + Parameters + ---------- + filename : str + Path to the .fits file + channel : int, default=0 + Channel of the image to access + """ + + def __init__(self, filename, channel=0): + self._fits_file = filename + self._channel = channel + From bf5e8e43db4d9079d375bf32a216b263271158c9 Mon Sep 17 00:00:00 2001 From: Jeff Jennings Date: Wed, 29 Nov 2023 17:32:29 -0500 Subject: [PATCH 3/7] ProcessFitsImage: add get_extent func --- src/mpol/input_output.py | 35 +++++++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) diff --git a/src/mpol/input_output.py b/src/mpol/input_output.py index 45093d4e..144cd26a 100644 --- a/src/mpol/input_output.py +++ b/src/mpol/input_output.py @@ -18,3 +18,38 @@ def __init__(self, filename, channel=0): self._fits_file = filename self._channel = channel + + def get_extent(self, header): + """Get extent (in RA and Dec, units of [arcsec]) of image""" + + # get the coordinate labels + nx = header["NAXIS1"] + ny = header["NAXIS2"] + + assert ( + nx % 2 == 0 and ny % 2 == 0 + ), f"Image dimensions x {nx} and y {ny} must be even." + + # RA coordinates + CDELT1 = 3600 * header["CDELT1"] # arcsec (converted from decimal deg) + # CRPIX1 = header["CRPIX1"] - 1.0 # Now indexed from 0 + + # DEC coordinates + CDELT2 = 3600 * header["CDELT2"] # arcsec + # CRPIX2 = header["CRPIX2"] - 1.0 # Now indexed from 0 + + RA = (np.arange(nx) - nx / 2) * CDELT1 # [arcsec] + DEC = (np.arange(ny) - ny / 2) * CDELT2 # [arcsec] + + # extent needs to include extra half-pixels. + # RA, DEC are pixel centers + + ext = ( + RA[0] - CDELT1 / 2, + RA[-1] + CDELT1 / 2, + DEC[0] - CDELT2 / 2, + DEC[-1] + CDELT2 / 2, + ) # [arcsec] + + return RA, DEC, ext + From c8de9764ba2dda4bd1f8f977377d9433ab2ef4bc Mon Sep 17 00:00:00 2001 From: Jeff Jennings Date: Wed, 29 Nov 2023 17:32:40 -0500 Subject: [PATCH 4/7] ProcessFitsImage: add get_beam func --- src/mpol/input_output.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/src/mpol/input_output.py b/src/mpol/input_output.py index 144cd26a..e5c3cf14 100644 --- a/src/mpol/input_output.py +++ b/src/mpol/input_output.py @@ -53,3 +53,23 @@ def get_extent(self, header): return RA, DEC, ext + + def get_beam(self, hdu_list, header): + """Get the major and minor widths [arcsec], and position angle, of a + clean beam""" + + if header.get("CASAMBM") is not None: + # Get the beam info from average of record array + data2 = hdu_list[1].data + BMAJ = np.median(data2["BMAJ"]) + BMIN = np.median(data2["BMIN"]) + BPA = np.median(data2["BPA"]) + else: + # Get the beam info from the header, like normal + BMAJ = 3600 * header["BMAJ"] + BMIN = 3600 * header["BMIN"] + BPA = header["BPA"] + + return BMAJ, BMIN, BPA + + From 48128ee9b359a62c78c7fca44118ed6cb1f6b250 Mon Sep 17 00:00:00 2001 From: Jeff Jennings Date: Wed, 29 Nov 2023 17:32:53 -0500 Subject: [PATCH 5/7] ProcessFitsImage: add get_image func --- src/mpol/input_output.py | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/src/mpol/input_output.py b/src/mpol/input_output.py index e5c3cf14..3ff1be3e 100644 --- a/src/mpol/input_output.py +++ b/src/mpol/input_output.py @@ -73,3 +73,29 @@ def get_beam(self, hdu_list, header): return BMAJ, BMIN, BPA + def get_image(self, beam=True): + """Load a .fits image and return as a numpy array. Also return image + extent and optionally (`beam`) the clean beam dimensions""" + + hdu_list = fits.open(self._fits_file) + hdu = hdu_list[0] + + if len(hdu.data.shape) in [3, 4]: + image = hdu.data[self._channel] # first channel + else: + image = hdu.data + + image *= 1e3 + + if len(image.shape) == 3: + image = np.squeeze(image) + + header = hdu.header + + RA, DEC, ext = self.get_extent(hdu.header) + + if beam: + return image, ext, self.get_beam(hdu_list, header) + else: + return image, ext + \ No newline at end of file From 21ea3d5cad75ae1daf87b6bd007c47a2b29b99f7 Mon Sep 17 00:00:00 2001 From: Jeff Jennings Date: Wed, 29 Nov 2023 20:18:17 -0500 Subject: [PATCH 6/7] add new test script --- test/input_output_test.py | 6 ++++++ 1 file changed, 6 insertions(+) create mode 100644 test/input_output_test.py diff --git a/test/input_output_test.py b/test/input_output_test.py new file mode 100644 index 00000000..e390b301 --- /dev/null +++ b/test/input_output_test.py @@ -0,0 +1,6 @@ +import numpy as np + +from astropy.utils.data import download_file + +from mpol.input_output import ProcessFitsImage + From 8a7340b455d70cd7b446f32f985da23ba3c89c84 Mon Sep 17 00:00:00 2001 From: Jeff Jennings Date: Wed, 29 Nov 2023 20:18:32 -0500 Subject: [PATCH 7/7] IO tests: add test_ProcessFitsImage --- test/input_output_test.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/test/input_output_test.py b/test/input_output_test.py index e390b301..39d1a4a6 100644 --- a/test/input_output_test.py +++ b/test/input_output_test.py @@ -4,3 +4,14 @@ from mpol.input_output import ProcessFitsImage +def test_ProcessFitsImage(): + # get a .fits file produced with casa + fname = download_file( + "https://zenodo.org/record/4711811/files/logo_cube.tclean.fits", + cache=True, + show_progress=True, + pkgname="mpol", + ) + + fits_image = ProcessFitsImage(fname) + clean_im, clean_im_ext, clean_beam = fits_image.get_image(beam=True) \ No newline at end of file