From 7eb0da79a54199817c476575a585761cc94d94cb Mon Sep 17 00:00:00 2001 From: Anno Knierim Date: Fri, 11 Oct 2024 14:14:16 +0200 Subject: [PATCH 1/6] Add analysis submodule with dynamic range and rms funtions --- radiotools/analysis/analysis.py | 74 +++++++++++++++++++++++++++++++++ 1 file changed, 74 insertions(+) create mode 100644 radiotools/analysis/analysis.py diff --git a/radiotools/analysis/analysis.py b/radiotools/analysis/analysis.py new file mode 100644 index 0000000..f61adcd --- /dev/null +++ b/radiotools/analysis/analysis.py @@ -0,0 +1,74 @@ +import torch +from astropy.io import fits +from numpy.typing import ArrayLike + +from radiotools.utils import rms + + +def get_source_rms(image: ArrayLike, center: tuple[int, int], *, offset=75): + """Get the root mean square (RMS) of a given source. + + Parameters + ---------- + image : array_like, shape (N, M) + Array containing image data of the source. + center : tuple[int, int] + Tuple of the center coordinates where the RMS + is calculated for a region spanned by ``offset``. + offset : int, optional + Offset that spans a region around ``center``: + + .. code:: + x0 = center_x - offset + x1 = center_x + offset + + y0 = center_y - offset + y1 = center_y + offset + + Default: 75 + + Returns + ------- + rms : float + RMS for the given source. + """ + center_x = int(center[0]) + center_y = int(center[1]) + + x0 = center_x - offset + x1 = center_x + offset + + y0 = center_y - offset + y1 = center_y + offset + + _rms = torch.from_numpy(rms(image[x0:x1, y0:y1])) + + return torch.mean(_rms) + + +def dynamic_range(path: str, *, offset: int = 75): + """Computes the dynamic range of a source. + + Parameters + ---------- + path : str or Path + Path to a FITS file containing image data of the source. + offset : int, optional + Offset for the region around the center pixel of the source. + Used to compute the RMS value of the source. Default: 75 + + Returns + ------- + dr : float + Dynamic range of the source. + """ + hdu = fits.open(path) + + image = hdu[0].data[0, 0, ...] + center = hdu[0].header["CRPIX1"], hdu[0].header["CRPIX2"] + + _rms = get_source_rms(image, center, offset=75) + + dr = image.max() / _rms + + return dr From eb86d07c4ff0ca60ec5eff6f245cbba0f0edb978 Mon Sep 17 00:00:00 2001 From: Anno Knierim Date: Fri, 11 Oct 2024 14:15:00 +0200 Subject: [PATCH 2/6] Add rms function to utils --- radiotools/utils/utils.py | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/radiotools/utils/utils.py b/radiotools/utils/utils.py index d6ffe63..229a4ac 100644 --- a/radiotools/utils/utils.py +++ b/radiotools/utils/utils.py @@ -1,5 +1,7 @@ +import numpy as np import requests from bs4 import BeautifulSoup +from numpy.typing import ArrayLike def get_array_names(url: str) -> list[str]: @@ -30,3 +32,24 @@ def get_array_names(url: str) -> list[str]: layouts = list(set(layouts)) return layouts + + +def rms(a: ArrayLike, *, axis: int | None = 0): + """Return an array of the root-mean-square (RMS) value of + the passed array. + + Parameters + ---------- + a : array_like + Array of which to find the RMS. + axis : int or None + + Returns + ------- + rms : np.ndarray + Array of rms values. + """ + if np.ndim(a) == 0: + axis = None + + return np.sqrt(np.mean(a**2, axis=axis)) From 0022862e0ba8363d4fa81796e79a65022812b011 Mon Sep 17 00:00:00 2001 From: Anno Knierim Date: Fri, 11 Oct 2024 14:15:30 +0200 Subject: [PATCH 3/6] Add new functions to inits --- radiotools/analysis/__init__.py | 3 +++ radiotools/utils/__init__.py | 4 ++-- 2 files changed, 5 insertions(+), 2 deletions(-) create mode 100644 radiotools/analysis/__init__.py diff --git a/radiotools/analysis/__init__.py b/radiotools/analysis/__init__.py new file mode 100644 index 0000000..d49666d --- /dev/null +++ b/radiotools/analysis/__init__.py @@ -0,0 +1,3 @@ +from .analysis import dynamic_range, get_source_rms + +__all__ = ["dynamic_range", "get_source_rms"] diff --git a/radiotools/utils/__init__.py b/radiotools/utils/__init__.py index 1720b0d..4197f3b 100644 --- a/radiotools/utils/__init__.py +++ b/radiotools/utils/__init__.py @@ -1,3 +1,3 @@ -from .utils import get_array_names +from .utils import get_array_names, rms -__all__ = ["get_array_names"] +__all__ = ["get_array_names", "rms"] From b66db5f7588ebe60809da044a9a2b14ea6b81162 Mon Sep 17 00:00:00 2001 From: Anno Knierim Date: Fri, 11 Oct 2024 14:15:52 +0200 Subject: [PATCH 4/6] Add tests for radiotools.analysis --- tests/test_analysis.py | 46 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) create mode 100644 tests/test_analysis.py diff --git a/tests/test_analysis.py b/tests/test_analysis.py new file mode 100644 index 0000000..54383cf --- /dev/null +++ b/tests/test_analysis.py @@ -0,0 +1,46 @@ +import numpy as np +from numpy.testing import assert_allclose + +from radiotools.utils import rms + + +class TestRMS: + def test_1d_array(self): + """Test 1d array.""" + a = np.array([10, 0, 150, 40, 30, 20, -50, -60, 5]) + expected = np.sqrt(1 / 9 * 31625) + + assert_allclose(rms(a), expected) + + def test_2d_array_axis_none(self): + """Test 2d array with axis=None.""" + a = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) + expected = 5.627314338711377 + + assert_allclose(rms(a, axis=None), expected) + + def test_2d_array_axis_0(self): + """Test 2d array with axis=0.""" + a = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) + expected = np.array( + [ + 4.69041575982343, + 5.5677643628300215, + 6.48074069840786, + ] + ) + + assert_allclose(rms(a, axis=0), expected) + + def test_2d_array_axis_1(self): + """Test 2d array with axis=1.""" + a = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) + expected = np.array( + [ + 2.160246899469287, + 5.066228051190222, + 8.04155872120988, + ] + ) + + assert_allclose(rms(a, axis=1), expected) From b63c870363638495ce39c25b046c7b2523948774 Mon Sep 17 00:00:00 2001 From: Anno Knierim Date: Fri, 11 Oct 2024 14:20:37 +0200 Subject: [PATCH 5/6] Add changelog --- docs/changes/8.feature.rst | 5 +++++ 1 file changed, 5 insertions(+) create mode 100644 docs/changes/8.feature.rst diff --git a/docs/changes/8.feature.rst b/docs/changes/8.feature.rst new file mode 100644 index 0000000..68bf6aa --- /dev/null +++ b/docs/changes/8.feature.rst @@ -0,0 +1,5 @@ +- Add analysis submodule + - Add ``radiotools.analysis.get_source_rms`` function + - Add ``radiotools.analysis.dynamic_range`` function that calculates the dynamic range as a function of the maximum pixel value of an image and the RMS of the source + +- Add ``rms`` function to ``radiotools.utils`` From aa02d39c2802f3210ca8fc5d95cfeaa4a9eaaf58 Mon Sep 17 00:00:00 2001 From: Anno Knierim Date: Thu, 14 Nov 2024 14:42:32 +0100 Subject: [PATCH 6/6] Test input with ndim=0 for radiotools.utils.rms --- tests/test_analysis.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/tests/test_analysis.py b/tests/test_analysis.py index 54383cf..46ec70e 100644 --- a/tests/test_analysis.py +++ b/tests/test_analysis.py @@ -44,3 +44,11 @@ def test_2d_array_axis_1(self): ) assert_allclose(rms(a, axis=1), expected) + + def test_input_ndim_0(self): + """Test input with ndim=0""" + a = np.int16(10) + + expected = 10.0 + + assert rms(a) == expected