Fix for units and numpy ops (#73)

* Fix for units and numpy ops

* Fix for readonly array

* New test

Author: AlecThomson

racs_tools/beamcon_2D.py

@@ -324,7 +324,7 @@ def get_common_beam(
     flags = np.array([beam == ZERO_BEAM for beam in beams]) | flags
 
     if cutoff is not None:
-        flags = beams.major > cutoff * u.arcsec | flags
+        flags = beams.major.to(u.arcsec).value > cutoff | flags
         if np.all(flags):
             logger.critical(
                 "All beams are larger than cutoff. All outputs will be blanked!"

racs_tools/convolve_uv.py

@@ -536,7 +536,7 @@ def smooth(
     if np.any(np.isnan(new_image)):
         logger.warning(f"{np.isnan(new_image).sum()} NaNs present in smoothed output")
 
-    new_image *= fac
+    new_image = fac * new_image
 
     # Ensure the output data-type is the same as the input
     return new_image.astype(out_dtype)

tests/test_2d.py

@@ -76,6 +76,55 @@ def make_2d_image(tmpdir) -> TestImage:
     outf.unlink()
 
 
+@pytest.fixture
+def make_2d_image_smaller(tmpdir) -> TestImage:
+    """Make a fake 2D image from with a Gaussian beam.
+
+    Args:
+        beam (Beam): Gaussian beam.
+
+    Returns:
+        str: Name of the output FITS file.
+    """
+    pix_scale = 2.5 * u.arcsec
+
+    beam = Beam(10 * u.arcsec, 10 * u.arcsec, 10 * u.deg)
+
+    data = beam.as_kernel(pixscale=pix_scale, x_size=100, y_size=100).array
+    data /= data.max()
+
+    hdu = fits.PrimaryHDU(data=data)
+
+    hdu.header = beam.attach_to_header(hdu.header)
+    hdu.header["BUNIT"] = "Jy/beam"
+    hdu.header["CDELT1"] = -pix_scale.to(u.deg).value
+    hdu.header["CDELT2"] = pix_scale.to(u.deg).value
+    hdu.header["CRPIX1"] = 50
+    hdu.header["CRPIX2"] = 50
+    hdu.header["CRVAL1"] = 0
+    hdu.header["CRVAL2"] = 0
+    hdu.header["CTYPE1"] = "RA---SIN"
+    hdu.header["CTYPE2"] = "DEC--SIN"
+    hdu.header["CUNIT1"] = "deg"
+    hdu.header["CUNIT2"] = "deg"
+    hdu.header["EQUINOX"] = 2000.0
+    hdu.header["RADESYS"] = "FK5"
+    hdu.header["LONPOLE"] = 180.0
+    hdu.header["LATPOLE"] = 0.0
+
+    outf = Path(tmpdir) / "2d_smaller.fits"
+    hdu.writeto(outf, overwrite=True)
+
+    yield TestImage(
+        path=outf,
+        beam=beam,
+        data=data,
+        pix_scale=pix_scale,
+    )
+
+    outf.unlink()
+
+
 @pytest.fixture
 def mirsmooth(make_2d_image: TestImage) -> TestImage:
     """Smooth a FITS image to a target beam using MIRIAD.
@@ -219,3 +268,14 @@ def test_smooth(make_2d_image: TestImage, mirsmooth: TestImage):
         assert np.allclose(
             smooth_data, mirsmooth.data, atol=1e-5
         ), f"Smooth with {conv_mode} does not match miriad"
+
+
+def test_get_common_beam(make_2d_image, make_2d_image_smaller):
+    common_beam = beamcon_2D.smooth_fits_files(
+        infile_list=[make_2d_image.path, make_2d_image_smaller.path],
+        cutoff=10,
+    )
+    assert (
+        common_beam.major.to(u.arcsec).value
+        == make_2d_image_smaller.beam.major.to(u.arcsec).value
+    )