PSF and Dirty image norm issue (or misuse)

[JulienNGirard]

Hello

With a colleague of mine, we are currently working on a new imager using part of codex-africanus code (which serves as an API for us, providing very useful methods for the gridding/degridding and other utility tasks).

• We are using Ducc's wgridder implementation (the wrapper to Ducc) after stop using nifty).
• We have been experiencing curious anomalies while inspecting the values from dirty generated in our code.
• As a matter of narrowing down the issue (or misuse...), I reproduced below a minimal example showing what is puzzling us.

On well-known real test datasets (from NenuFAR) for which we have a good knowledge of, we tried to generate dirty images (and PSF following the same track as Simon due to the lack of a dedicated PSF method.

Here is the code we have used to generate the dirty

from daskms import xds_from_ms,xds_from_table
from africanus.gridding.wgridder import dirty as dirty_np

# loading MS
data_path = data_folder / 'SB155.rebin.MS'
data_path = data_path.resolve()
ms_file = str(data_path)
ms_file_ant = str(data_path / "ANTENNA")
ms = xds_from_ms(ms_file)[0]

# getting required information from MS

ms_spectral_window = xds_from_table(ms_file + "/SPECTRAL_WINDOW")[0]
freq = ms_spectral_window.REF_FREQUENCY.compute().to_numpy()                  # .values() is shorter... sorry
freq_bin_idx = np.array([0])
freq_bin_counts = np.array([1])

uvw = ms.UVW.data.compute()
vis = (ms.DATA.isel(corr=0).data+ms.DATA.isel(corr=-1).data).compute()         # Stokes I from DATA column

# image parameters
nx = 1024
ny = 1024
wgt = ms.WEIGHT_SPECTRUM.isel(corr=0).data.compute()                                # taking from XX (same from YY)
cell= np.radians(5./60)

nthreads = 1

# creating dirty 
image = dirty_np(
    uvw,
    freq,
    vis,
    freq_bin_idx,
    freq_bin_counts,
    nx,
    ny,
    cell=cell,
    epsilon=1e-06,
    weights=wgt,
    nthreads=nthreads,
    do_wstacking=False
)

With this code, we tried 3 different cases.

(the Flag options does not seem to be working for us, but it will be another issue as we do not understand what is the returned error).

1. Virgo A dirty:

Dirty Image produced with the Codex Africanus code (forget about the missing origin="lower" in the plot)

Dirty image produced by WSClean

Corresponding Wsclean command (set to natural weighting):

wsclean -make-psf -mem 95 -no-reorder -no-update-model-required -mgain 0.8 -weight briggs 2 -size 1024 1024 -scale 5arcmin -pol I -data-column DATA -niter 0 -name testdirtySB155-2 SB155.rebin.MS

Comparison:

• Codex Africanus Dirty peak flux is ~1e15 Jy/beam
• WSclean Dirty peak flux is ~1e7 Jy/beam

On this dataset, we have been fighting with this 1e7 factor ratio in flux without being able to understand what is happening.
Is that anything to do with some kind of normalization that should be taken into account? corrected for?

2. PSF generation:

PSF also do not share the same norm.

• WSClean properly normalized to 1.
• Codex: using unity visibilities and the (natural) weighting (see here)
  Question:The normalized PSF is also high in power (1e5 - 1e8) on some dataset but is not consistent with the ratio we measure on dirty images between WSclean and Codex Africanus.

3. Cygnus loop (Dirty) with NenuFAR:

• WSClean dirty of CORRECTED_DATA: Max = 119 Jy/beam

- Codex Africanus dirty of CORRECTED_DATA: Max = 9.7e13 Jy/beam

Ratio=1e12 for this dataset.

Few remarks:

• we tried to only keep the unflagged data by filtering them manually out of the vis table and feeding them to dirty_np() did not solve the issue, so it is not due to bad data.
• we have tried to change he image pixel size, but it did not change the ratio value for a given dataset
• the ratio seems to vary from one dataset to another
• tried squareroot(weights) & square(weights) to see if it had any impact, the ratio value changed a bit but not in a large extent.

Is there some normalisation we are missing? weights? Nvis? FFT norm? Amount for flagged/unflagged data?

And also in the documentation (about the dirty method: here what is meant by "prewhitening" the visibility data? Is that what is missing?

Please help, we tried many things without being able to pinpoint what is wrong in our usage of the method.

Many thanks

[landmanbester]

Hi @JulienNGirard. These wrappers might be slightly outdated, that's all on me, sorry. I've been working on some stuff that might be relevant to you but hasn't quite made it into codex-africanus. Do you have time for a video call? It would be great to catch up

[landmanbester]

Ah, I see the issue. There is no normalisation by wsum in these wrappers, you have to take care of that yourself. If you are not using any weights then wsum is just the number of data points going into the gridder otherwise it is the sum of the weights. The peak of the PSF should be equal to this value but due to gridding precision it may be slight off (that's the reason for the double_accum option). Not sure if you are aware of https://github.com/ratt-ru/pfb-imaging. There are some gridding utilities that might be useful to you (see e.g. https://github.com/ratt-ru/pfb-imaging/blob/main/pfb/operators/gridder.py and https://github.com/ratt-ru/pfb-imaging/blob/main/pfb/utils/weighting.py)

[JulienNGirard]

Thank very much Landman! (long time no see, going off for holiday tonight!) :)) thank you for the feedback, I suspected that it was some normalisation issue.
We were using the "IMAGING_WEIGHT" column as input.
We are going to test that right away.

We are also considering using pfb-imaging as an alternate plan, thanks to your input !