"""
Functions that define and manipulate ConvolutionFunctions.
The griddata has axes [chan, pol, z, dy, dx, y, x] where z, y, x
are spatial axes in either sky or Fourier plane. The order in the
WCS is reversed so the grid_WCS describes UU, VV, DUU, DVV, WW, STOKES, FREQ axes.
GridData can be used to hold the Fourier transform of an Image or gridded
visibilities. In addition, the convolution function can be stored in a
GridData, most probably with finer spatial sampling.
"""
__all__ = [
"calculate_bounding_box_convolutionfunction",
"apply_bounding_box_convolutionfunction",
"export_convolutionfunction_to_fits",
]
import logging
import numpy
from astropy.io import fits
from ska_sdp_datamodels.gridded_visibility.grid_vis_model import ConvolutionFunction
log = logging.getLogger("rascil-logger")
[docs]
def apply_bounding_box_convolutionfunction(cf, fractional_level=1e-4):
"""Apply a bounding box to a convolution function
:param cf:
:param fractional_level:
:return: bounded convolution function
"""
newcf = cf.copy(deep=True)
nx = newcf["pixels"].data.shape[-1]
ny = newcf["pixels"].data.shape[-2]
mask = numpy.max(numpy.abs(newcf["pixels"].data), axis=(0, 1, 2, 3, 4))
mask /= numpy.max(mask)
coords = numpy.argwhere(mask > fractional_level)
crpx = int(numpy.round(cf.convolutionfunction_acc.cf_wcs.wcs.crpix[0]))
crpy = int(numpy.round(cf.convolutionfunction_acc.cf_wcs.wcs.crpix[1]))
x0, y0 = coords.min(axis=0, initial=cf["pixels"].data.shape[-1])
dx = crpx - x0
dy = crpy - y0
x0 -= 1
y0 -= 1
x1 = crpx + dx - 1
y1 = crpy + dy - 1
nny, nnx = newcf["pixels"].data.shape[-2], newcf["pixels"].data.shape[-1]
newcf.convolutionfunction_acc.cf_wcs.wcs.crpix[0] += nnx / 2 - nx / 2
newcf.convolutionfunction_acc.cf_wcs.wcs.crpix[1] += nny / 2 - ny / 2
newcf = ConvolutionFunction.constructor(
data=newcf["pixels"].data[..., y0:y1, x0:x1],
cf_wcs=newcf.convolutionfunction_acc.cf_wcs,
polarisation_frame=newcf.convolutionfunction_acc.polarisation_frame,
)
return newcf
[docs]
def calculate_bounding_box_convolutionfunction(cf, fractional_level=1e-4):
"""Calculate bounding boxes
Returns a list of bounding boxes where each element is
(z, (y0, y1), (x0, x1))
These can be used in griddata/degridding.
:param cf:
:param fractional_level:
:return: list of bounding boxes
"""
bboxes = list()
threshold = fractional_level * numpy.max(numpy.abs(cf["pixels"].data))
for z in range(cf["pixels"].data.shape[2]):
mask = numpy.max(numpy.abs(cf["pixels"].data[:, :, z, ...]), axis=(0, 1, 2, 3))
coords = numpy.argwhere(mask > threshold)
x0, y0 = coords.min(axis=0, initial=cf["pixels"].data.shape[-1])
x1, y1 = coords.max(axis=0, initial=cf["pixels"].data.shape[-1])
bboxes.append((z, (y0, y1), (x0, x1)))
return bboxes
[docs]
def export_convolutionfunction_to_fits(
cf: ConvolutionFunction, fitsfile: str = "cf.fits"
):
"""Write a convolution function to fits
:param cf: Convolu
:param fitsfile: Name of output fits file in storage
:returns: None
See also
:py:func:`rascil.processing_components.image.operations.import_image_from_array`
"""
return fits.writeto(
filename=fitsfile,
data=numpy.real(cf["pixels"].data),
header=cf.convolutionfunction_acc.cf_wcs.to_header(),
overwrite=True,
)