diff --git a/CHANGELOG.md b/CHANGELOG.md index d3a12595..247fb9a9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -29,6 +29,8 @@ - removing the `pol` string in the polarisation field of the `ProcessedNameComponents` - `wclean` output `-` separation character chhanged to `.` +- The median RMS of the field image is calculated on an inner region when + cnnstructing the validation plot # 0.2.5 diff --git a/flint/coadd/linmos.py b/flint/coadd/linmos.py index 307466f7..50ac2199 100644 --- a/flint/coadd/linmos.py +++ b/flint/coadd/linmos.py @@ -41,7 +41,8 @@ class BoundingBox(NamedTuple): def create_bound_box(image_data: np.ndarray, is_masked: bool = False) -> BoundingBox: - """Construct a bounding box around finite pixels. + """Construct a bounding box around finite pixels for a 2D image. This does not + support cube type images. If ``is_mask` is ``False``, the ``image_data`` will be masked internally using ``numpy.isfinite``. diff --git a/flint/validation.py b/flint/validation.py index a85a1510..42acc074 100644 --- a/flint/validation.py +++ b/flint/validation.py @@ -203,6 +203,37 @@ class MatchResult(NamedTuple): """Brightness in Jy of source in the second survey""" +def extract_inner_image_array_region(image: np.ndarray, fraction: float) -> np.ndarray: + """Extract an inner region of an image array. The size of the extracted array + is specified via the `fraction` parameter, which should be in the range of 0 to 1. + + Args: + image (np.ndarray): The image to extract pixels from + fraction (float): The size of the region to extraction. + + Returns: + np.ndarray: The extracted sub-region + """ + assert 0.0 < fraction <= 1.0, f"{fraction=}, but has to be between 0.0 to 1.0" + + image_shape = image.shape + yx_shape = np.array(image_shape[-2:]) + + yx_center = np.ceil(yx_shape / 2).astype(int) + box_size = (yx_shape * fraction).astype(int) + box_delta = np.ceil(box_size / 2).astype(int) + + min_y = max(0, yx_center[0] - box_delta[0]) + max_y = min(yx_shape[0], yx_center[0] + box_delta[0]) + + min_x = max(0, yx_center[1] - box_delta[1]) + max_x = min(yx_shape[1], yx_center[1] + box_delta[1]) + + sub_image_array = image[..., min_y:max_y, min_x:max_x] + + return sub_image_array + + def get_known_catalogue_info(name: str) -> Catalogue: """Return the parameters of a recognised catalogue. @@ -291,11 +322,16 @@ def load_known_catalogue( return table, catalogue -def get_rms_image_info(rms_path: Path) -> RMSImageInfo: +def get_rms_image_info(rms_path: Path, extract_fraction: float = 0.5) -> RMSImageInfo: """Extract information about the RMS image and construct a representative structure + When computing the pixel statistics an inner region of the RMS map is extracted + to avoid the primary beam roll off. Other properties of the ``RMSImageInfo`` are + computed on the whole image (e.g. area) + Args: rms_path (Path): The RMS image that will be presented + extract_fraction (float, optional): Extract the inner region of the base RMS image map. Defaults to 0.5. Returns: RMSImageInfo: Extracted RMS image information @@ -305,6 +341,11 @@ def get_rms_image_info(rms_path: Path) -> RMSImageInfo: rms_header = rms_image[0].header # type: ignore rms_data = rms_image[0].data # type: ignore + logger.info(f"Extracting inner region of {extract_fraction=} for RMS statistics") + sub_rms_data = extract_inner_image_array_region( + image=rms_data, fraction=extract_fraction + ) + valid_pixels = np.sum(np.isfinite(rms_data) & (rms_data != 0)) ra_cell_size = rms_header["CDELT1"] @@ -324,11 +365,11 @@ def get_rms_image_info(rms_path: Path) -> RMSImageInfo: area=area, wcs=wcs, centre=centre_sky, - median=np.nanmedian(rms_data), - minimum=np.nanmin(rms_data), - maximum=np.nanmax(rms_data), - std=float(np.nanstd(rms_data)), - mode=stats.mode(rms_data, keepdims=False).mode[0], + median=float(np.nanmedian(sub_rms_data)), + minimum=float(np.nanmin(sub_rms_data)), + maximum=float(np.nanmax(sub_rms_data)), + std=float(np.nanstd(sub_rms_data)), + mode=stats.mode(sub_rms_data, keepdims=False).mode[0], ) return rms_image_info diff --git a/tests/test_linmos_coadd.py b/tests/test_linmos_coadd.py index fd4b8d6f..a3968a0d 100644 --- a/tests/test_linmos_coadd.py +++ b/tests/test_linmos_coadd.py @@ -119,6 +119,7 @@ def test_trim_fits_image_matching(tmp_path): def test_bounding_box(): + """Create a bounding box around a region of valid non-nan/inf pixels""" data = np.zeros((1000, 1000)) data[10:600, 20:500] = 1 data[data == 0] = np.nan @@ -132,7 +133,18 @@ def test_bounding_box(): assert bb.ymax == 499 # slices upper limit no inclusive +def test_bounding_box_cube(): + """Cube cut bounding boxes. Currently not supported.""" + data = np.zeros((3, 1000, 1000)) + data[:, 10:600, 20:500] = 1 + data[data == 0] = np.nan + + with pytest.raises(AssertionError): + create_bound_box(image_data=data) + + def test_bounding_box_with_mask(): + """Create a bounding box where the input is converted to a boolean arrau""" data = np.zeros((1000, 1000)) data[10:600, 20:500] = 1 diff --git a/tests/test_validation.py b/tests/test_validation.py index 9a755e77..eb151053 100644 --- a/tests/test_validation.py +++ b/tests/test_validation.py @@ -13,6 +13,7 @@ RMSImageInfo, SourceCounts, calculate_area_correction_per_flux, + extract_inner_image_array_region, get_parser, get_rms_image_info, get_source_counts, @@ -48,6 +49,24 @@ def comp_path(tmpdir): return comp_path +def test_extract_inner_image_array_region(): + """Get the inner region of an image array""" + shape = (100, 100) + image = np.arange(np.prod(shape)).reshape(shape) + + sub_image = extract_inner_image_array_region(image=image, fraction=0.5) + assert sub_image.shape == (50, 50) + + shape = (1, 1, 100, 100) + image = np.arange(np.prod(shape)).reshape(shape) + + sub_image = extract_inner_image_array_region(image=image, fraction=0.5) + assert sub_image.shape == (1, 1, 50, 50) + + with pytest.raises(AssertionError): + extract_inner_image_array_region(image=image, fraction=22.5) + + def test_plot_source_counts(comp_path, rms_path): """Just a simple test to make sure there are no immediate matplotlib api errors. Not testing for actual plotting results.""" @@ -96,14 +115,16 @@ def test_calculate_area_correction(rms_path): def test_rms_image_info(rms_path): rms_info = get_rms_image_info(rms_path=rms_path) + print(rms_info) + assert isinstance(rms_info, RMSImageInfo) assert rms_info.path == rms_path assert rms_info.no_valid_pixels == 150 assert rms_info.shape == (10, 15) - assert np.isclose(0.0001515522, rms_info.median) - assert np.isclose(0.00015135764, rms_info.minimum) - assert np.isclose(0.0001518184, rms_info.maximum) - assert np.isclose(1.1098655e-07, rms_info.std) + assert np.isclose(0.0001515522, rms_info.median, atol=1e-5) + assert np.isclose(0.00015135764, rms_info.minimum, atol=1e-5) + assert np.isclose(0.0001518184, rms_info.maximum, atol=1e-5) + assert np.isclose(1.1098655e-07, rms_info.std, atol=1e-5) class Example(NamedTuple):