This is an enhancement rather than a fix, but to play better with other code I moved the print() commands in beacon_2D to use astropy.log and set the verbosity level appropriately using the command-line arguments. The reason I'm mentioning this is because my fork:
https://github.com/dlakaplan/RACS-tools
has those of those fixes alongside the other two issues I raised, and I wanted to be clear what the motivation was. I personally find it easier this way, but others may have different opinions.

In any case, if you want to merge in this or the other issues I raised (which are more in the ways of fixes) let me know and I'll open a PR. I can probably try to separate the different issues if I need to.

Thanks for the useful tool. When trying to use it I found that the getimdata routine assumes that all images are 4D (including polarization and spectral axes). So I modified it in my fork:
https://github.com/dlakaplan/RACS-tools
to work on 2D (squeezed) images as well.

As a related issue, the determination of the nx,ny dimensions was backwards as far as I understand (although for square images it doesn't matter): the images are stored as [NY,NX]. So I changed that too.

I was doing this with @marxide

Please make!

Reported by @axshen - mosaic images can contain NaNs. Especially on the border, but potentially in cases of missing beams. We should work to make our robust method more robust to NaNs.

We need to update conolve_uv.py.

Our approach will be to convert NaN to 0, fix boundary pixels, and reapply NaNs in the output image. @wasimraja81 has a method using iterative clipping in the Fourier domain to find and fix edge pixels that will be affected.

Testing out getnoise_list.py on a MeerKAT polarisation cube. I've constructed Q and U frequency cubes using WSclean which I've convolved to a common resolution and stacked into a RA / Dec / freq cube. On trying to run getnoise_list I encounter a UnitConversionError:

$ python3 getnoise_list.py --file bad_channels_chan512.txt --ncores 32 RMSynth_chan512/qcube.fits RMSynth_chan512/ucube.fits
Checking Q: 100%|█████████████████████████████████| 443/443 [00:01<00:00, 286.68it/s]
Traceback (most recent call last):
  File "/homes/riseley/.local/lib/python3.9/site-packages/astropy/units/quantity.py", line 890, in to_value
    scale = self.unit._to(unit)
  File "/homes/riseley/.local/lib/python3.9/site-packages/astropy/units/core.py", line 1093, in _to
    raise UnitConversionError(
astropy.units.core.UnitConversionError: 'Unit("Jy / beam")' is not a scaled version of 'Unit(dimensionless)'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/homes/riseley/.local/lib/python3.9/site-packages/astropy/units/quantity.py", line 1248, in __float__
    return float(self.to_value(dimensionless_unscaled))
  File "/homes/riseley/.local/lib/python3.9/site-packages/astropy/units/quantity.py", line 893, in to_value
    value = self._to_value(unit, equivalencies)
  File "/homes/riseley/.local/lib/python3.9/site-packages/astropy/units/quantity.py", line 800, in _to_value
    return self.unit.to(unit, self.view(np.ndarray),
  File "/homes/riseley/.local/lib/python3.9/site-packages/astropy/units/core.py", line 1130, in to
    return self._get_converter(Unit(other),
  File "/homes/riseley/.local/lib/python3.9/site-packages/astropy/units/core.py", line 1061, in _get_converter
    raise exc
  File "/homes/riseley/.local/lib/python3.9/site-packages/astropy/units/core.py", line 1046, in _get_converter
    return self._apply_equivalencies(
  File "/homes/riseley/.local/lib/python3.9/site-packages/astropy/units/core.py", line 1022, in _apply_equivalencies
    raise UnitConversionError(
astropy.units.core.UnitConversionError: 'Jy / beam' and '' (dimensionless) are not convertible

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/local/work/riseley/software/RACS-tools/racs_tools/getnoise_list.py", line 330, in <module>
    cli()
  File "/local/work/riseley/software/RACS-tools/racs_tools/getnoise_list.py", line 324, in cli
    main(pool, args)
  File "/local/work/riseley/software/RACS-tools/racs_tools/getnoise_list.py", line 214, in main
    totalbad = getbadchans(
  File "/local/work/riseley/software/RACS-tools/racs_tools/getnoise_list.py", line 123, in getbadchans
    qnoisevals = np.array(qnoisevals)
  File "/homes/riseley/.local/lib/python3.9/site-packages/astropy/units/quantity.py", line 1250, in __float__
    raise TypeError('only dimensionless scalar quantities can be '
TypeError: only dimensionless scalar quantities can be converted to Python scalars

As a quick sanity check I removed the Jy/beam BUNIT from the Q/U-cube headers but that immediately threw an error when I tried running getnoise_list so it seems to be a required keyword.

Suggestion from Paddy Leahy - it's possible in some other convolvers to choose a resolution higher than the lowest common resolution. In this case, rather than convolving with a Gaussian, the kernal is an inverse Gaussian - upweighting longer baselines. The tradeoff is that the noise becomes larger. This can be acceptable within some reasonable range.

@wasimraja81 - we should discuss to see if this is possible in our analytic convolution scheme. Care would need to be taken to preserve units in Jy/beam.

Major flux scaling issues occur when the convolving beam is very small and therefore undersampled. This occurs when the target beam is very close to the original beam.

A similar issue was found in QUOCKA. Solution there was to calculate the critically sampled beam and up the target to the critical value.

Miriad deals with this somehow. Further investigation needed. A warning should at least be generated.

Another issue: while some images store the pixel scales via CDELT, we have some where the WCS has been replaced by the CD matrix. So I replaced the calls to get the pixel scales with a more general astropy routine:

        w = astropy.wcs.WCS(hdu[0])
        pixelscales = astropy.wcs.utils.proj_plane_pixel_scales(w)
        
        dxas = pixelscales[0]*u.deg
        dyas = pixelscales[1]*u.deg

in getimdata and similar in getmaxbeam.

Hi Alec:

I'm currently thinking about how to set up a convolution step for the VLASS-polarization pipeline, and this seems like it might be the least fuss to implement (the main alternative I'm considering being CASA). I used this tool with great success while prototyping some VLASS results last week.

My biggest complaint at the moment is that it's not immediately clear how easy it would be to integrate this into a Python pipeline (i.e., not using the command line interface). The smooth() function looks appropriate, but it's not clear what the inputs are. Would you consider making it more script-friendly (or at a minimum, commenting it a bit more thoroughly so I can see how to use the existing functions)? I'm mostly interested in the 2D version, since the VLASS data is packaged as single planes. In the ideal case (if you're feeling generous), I'm envisioning a function that can I feed an input array, input beam size, target beam size, and method; and get back a convolved array.

Thanks,
Cameron