Merge pull request #166 from tjgalvin/inneregion

Calculate RMS of inner region during validation plot

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
 

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``.
 

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

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
 

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):