Skip to content

Commit b2bfcee

Browse files
timcornwelltimcornwell
committed
Use convert iso coalesce
1 parent 75aefcf commit b2bfcee

19 files changed

Lines changed: 345 additions & 237 deletions

File tree

‎processing_components/calibration/calibration.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -23,7 +23,7 @@
2323

2424
from ..visibility.base import create_visibility_from_rows
2525
from ..calibration.operations import apply_gaintable, create_gaintable_from_blockvisibility
26-
from ..visibility.coalesce import convert_blockvisibility_to_visibility, decoalesce_visibility
26+
from ..visibility.coalesce import convert_blockvisibility_to_visibility, convert_visibility_to_blockvisibility
2727
from ..visibility.base import copy_visibility
2828
from ..imaging.base import predict_skycomponent_visibility, predict_2d
2929
from ..visibility.operations import divide_visibility

‎processing_components/imaging/base.py‎

Lines changed: 4 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -43,7 +43,7 @@
4343
degrid_visibility_from_griddata
4444
from ..griddata.operations import create_griddata_from_image
4545
from ..visibility.base import copy_visibility, phaserotate_visibility
46-
from ..visibility.coalesce import coalesce_visibility, decoalesce_visibility, convert_blockvisibility_to_visibility
46+
from ..visibility.coalesce import convert_blockvisibility_to_visibility, convert_visibility_to_blockvisibility
4747

4848
log = logging.getLogger(__name__)
4949

@@ -116,8 +116,7 @@ def predict_2d(vis: Union[BlockVisibility, Visibility], model: Image, gcfcf=None
116116
:return: resulting visibility (in place works)
117117
"""
118118
if isinstance(vis, BlockVisibility):
119-
log.debug("imaging.predict: coalescing prior to prediction")
120-
avis = coalesce_visibility(vis, **kwargs)
119+
avis = convert_blockvisibility_to_visibility(vis)
121120
else:
122121
avis = vis
123122

@@ -141,7 +140,7 @@ def predict_2d(vis: Union[BlockVisibility, Visibility], model: Image, gcfcf=None
141140

142141
if isinstance(vis, BlockVisibility) and isinstance(svis, Visibility):
143142
log.debug("imaging.predict decoalescing post prediction")
144-
return decoalesce_visibility(svis)
143+
return convert_visibility_to_blockvisibility(svis)
145144
else:
146145
return svis
147146

@@ -164,7 +163,7 @@ def invert_2d(vis: Visibility, im: Image, dopsf: bool = False, normalize: bool =
164163
165164
"""
166165
if not isinstance(vis, Visibility):
167-
svis = coalesce_visibility(vis, **kwargs)
166+
svis = convert_blockvisibility_to_visibility(vis)
168167
else:
169168
svis = copy_visibility(vis)
170169

‎processing_components/imaging/timeslice_single.py‎

Lines changed: 17 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -23,15 +23,10 @@
2323
2424
2525
"""
26-
import numpy
27-
from scipy.interpolate import RegularGridInterpolator
28-
from data_models.memory_data_models import Visibility, Image
29-
from ..image.operations import copy_image, create_empty_image_like
30-
from ..imaging.base import predict_2d, invert_2d
31-
from ..visibility.coalesce import coalesce_visibility
32-
3326
import logging
3427

28+
from ..visibility.coalesce import convert_blockvisibility_to_visibility, convert_visibility_to_blockvisibility
29+
3530
log = logging.getLogger(__name__)
3631

3732
"""
@@ -60,15 +55,14 @@
6055
6156
"""
6257
import numpy
63-
from scipy.interpolate import griddata, RectBivariateSpline, RegularGridInterpolator, interp2d
6458

65-
from data_models.memory_data_models import Visibility, Image
59+
from data_models.memory_data_models import BlockVisibility, Visibility, Image
6660

67-
from processing_components.image.operations import copy_image, create_empty_image_like, reproject_image
61+
from processing_components.image.operations import reproject_image
6862

6963
from processing_components.imaging.base import predict_2d, invert_2d
7064

71-
from processing_components.visibility.coalesce import coalesce_visibility
65+
from processing_components.visibility.coalesce import convert_blockvisibility_to_visibility
7266

7367

7468
def fit_uvwplane_only(vis: Visibility) -> (float, float):
@@ -124,10 +118,12 @@ def predict_timeslice_single(vis: Visibility, model: Image, predict=predict_2d,
124118
log.debug("predict_timeslice: predicting using time slices")
125119

126120
if not isinstance(vis, Visibility):
127-
avis = coalesce_visibility(vis, **kwargs)
121+
avis = convert_blockvisibility_to_visibility(vis)
128122
else:
129123
avis = vis
130-
124+
125+
avis.data['vis'][...] = 0.0
126+
131127
# Fit and remove best fitting plane for this slice
132128
uvw = avis.uvw
133129
avis, p, q = fit_uvwplane(avis, remove=remove)
@@ -145,10 +141,13 @@ def predict_timeslice_single(vis: Visibility, model: Image, predict=predict_2d,
145141
# Now we can do the predict
146142
if remove:
147143
avis.data['uvw'][...] = uvw
148-
144+
149145
avis = predict(avis, workimage, gcfcf=gcfcf, **kwargs)
150146

151-
return avis
147+
if isinstance(vis, BlockVisibility) and isinstance(avis, Visibility):
148+
return convert_visibility_to_blockvisibility(avis)
149+
else:
150+
return avis
152151

153152

154153
def invert_timeslice_single(vis: Visibility, im: Image, dopsf, normalize=True, remove=True,
@@ -163,7 +162,7 @@ def invert_timeslice_single(vis: Visibility, im: Image, dopsf, normalize=True, r
163162
:param normalize: Normalize by the sum of weights (True)
164163
"""
165164
if not isinstance(vis, Visibility):
166-
avis = coalesce_visibility(vis, **kwargs)
165+
avis = convert_blockvisibility_to_visibility(vis)
167166
else:
168167
avis = vis
169168

@@ -185,5 +184,5 @@ def invert_timeslice_single(vis: Visibility, im: Image, dopsf, normalize=True, r
185184

186185
if remove:
187186
avis.data['uvw'][...] = uvw
188-
189-
return finalimage, sumwt
187+
188+
return finalimage, sumwt

‎processing_components/imaging/wstack_single.py‎

Lines changed: 6 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -17,7 +17,7 @@
1717

1818
from ..image.operations import copy_image
1919
from ..visibility.base import copy_visibility
20-
from ..visibility.coalesce import coalesce_visibility, decoalesce_visibility
20+
from ..visibility.coalesce import convert_blockvisibility_to_visibility, convert_visibility_to_blockvisibility
2121
from ..imaging.base import predict_2d, invert_2d
2222

2323
import logging
@@ -36,10 +36,12 @@ def predict_wstack_single(vis, model, remove=True, gcfcf=None, **kwargs) -> Visi
3636

3737
if not isinstance(vis, Visibility):
3838
log.debug("predict_wstack_single: Coalescing")
39-
avis = coalesce_visibility(vis, **kwargs)
39+
avis = convert_blockvisibility_to_visibility(vis)
4040
else:
4141
avis = vis
4242

43+
avis.data['vis'][...] = 0.0
44+
4345
log.debug("predict_wstack_single: predicting using single w slice")
4446

4547
# We might want to do wprojection so we remove the average w
@@ -50,7 +52,7 @@ def predict_wstack_single(vis, model, remove=True, gcfcf=None, **kwargs) -> Visi
5052

5153
# Calculate w beam and apply to the model. The imaginary part is not needed
5254
workimage = copy_image(model)
53-
w_beam = create_w_term_like(model, w_average, vis.phasecentre)
55+
w_beam = create_w_term_like(model, w_average, avis.phasecentre)
5456

5557
# Do the real part
5658
workimage.data = w_beam.data.real * model.data
@@ -65,8 +67,7 @@ def predict_wstack_single(vis, model, remove=True, gcfcf=None, **kwargs) -> Visi
6567
avis.data['uvw'][..., 2] += w_average
6668

6769
if isinstance(vis, BlockVisibility) and isinstance(avis, Visibility):
68-
log.debug("imaging.predict decoalescing post prediction")
69-
return decoalesce_visibility(avis)
70+
return convert_visibility_to_blockvisibility(avis)
7071
else:
7172
return avis
7273

‎processing_components/simulation/testing_support.py‎

Lines changed: 91 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -40,22 +40,22 @@
4040
from astropy.wcs.utils import pixel_to_skycoord
4141
from scipy import interpolate
4242

43-
from data_models.memory_data_models import Configuration, Image, GainTable, Skycomponent
43+
from data_models.memory_data_models import Configuration, Image, GainTable, Skycomponent, SkyModel
4444
from data_models.parameters import arl_path
4545
from data_models.parameters import get_parameter
4646
from data_models.polarisation import PolarisationFrame
47-
from processing_library.image.operations import create_image_from_array
48-
from processing_library.util.coordinate_support import xyz_at_latitude
4947
from processing_components.calibration.calibration_control import create_calibration_controls
5048
from processing_components.calibration.operations import create_gaintable_from_blockvisibility, apply_gaintable
51-
from processing_components.image.operations import import_image_from_fits, qa_image
49+
from processing_components.image.operations import import_image_from_fits
5250
from processing_components.imaging.base import predict_2d, predict_skycomponent_visibility, \
5351
create_image_from_visibility, advise_wide_field
52+
from processing_components.imaging.primary_beams import create_pb
5453
from processing_components.skycomponent.operations import create_skycomponent, insert_skycomponent, \
5554
apply_beam_to_skycomponent, filter_skycomponents_by_flux
5655
from processing_components.visibility.base import create_blockvisibility, create_visibility
5756
from processing_components.visibility.coalesce import convert_visibility_to_blockvisibility
58-
from processing_components.imaging.primary_beams import create_pb
57+
from processing_library.image.operations import create_image_from_array
58+
from processing_library.util.coordinate_support import xyz_at_latitude
5959

6060
log = logging.getLogger(__name__)
6161

@@ -306,7 +306,7 @@ def create_test_image_from_s3(npixel=16384, polarisation_frame=PolarisationFrame
306306
assert fov in [10, 20, 40], "Field of view invalid: use one of %s" % ([10, 20, 40])
307307
csvfilename = arl_path('data/models/S3_151MHz_%ddeg.csv' % (fov))
308308
log.info('create_test_image_from_s3: Reading S3 sources from %s ' % csvfilename)
309-
309+
310310
with open(csvfilename) as csvfile:
311311
readCSV = csv.reader(csvfile, delimiter=',')
312312
r = 0
@@ -328,11 +328,11 @@ def create_test_image_from_s3(npixel=16384, polarisation_frame=PolarisationFrame
328328
r += 1
329329

330330
csvfile.close()
331-
331+
332332
assert len(fluxes) > 0, "No sources found above flux limit %s" % flux_limit
333-
333+
334334
log.info('create_test_image_from_s3: %d sources read' % (len(fluxes)))
335-
335+
336336
p = w.sub(2).wcs_world2pix(numpy.array(ras), numpy.array(decs), 1)
337337
total_flux = numpy.sum(fluxes)
338338
fluxes = numpy.array(fluxes)
@@ -361,9 +361,9 @@ def create_low_test_image_from_gleam(npixel=512, polarisation_frame=Polarisation
361361
"""Create LOW test image from the GLEAM survey
362362
363363
Stokes I is estimated from a cubic spline fit to the measured fluxes. The polarised flux is always zero.
364-
364+
365365
See http://www.mwatelescope.org/science/gleam-survey The catalog is available from Vizier.
366-
366+
367367
VIII/100 GaLactic and Extragalactic All-sky MWA survey (Hurley-Walker+, 2016)
368368
369369
GaLactic and Extragalactic All-sky Murchison Wide Field Array (GLEAM) survey. I: A low-frequency extragalactic
@@ -377,7 +377,7 @@ def create_low_test_image_from_gleam(npixel=512, polarisation_frame=Polarisation
377377
:param phasecentre: phasecentre (SkyCoord)
378378
:param kind: Kind of interpolation (see scipy.interpolate.interp1d) Default: linear
379379
:return: Image
380-
380+
381381
"""
382382

383383
if phasecentre is None:
@@ -418,6 +418,81 @@ def create_low_test_image_from_gleam(npixel=512, polarisation_frame=Polarisation
418418
return model
419419

420420

421+
def create_low_test_skymodel_from_gleam(npixel=512, polarisation_frame=PolarisationFrame("stokesI"), cellsize=0.000015,
422+
frequency=numpy.array([1e8]), channel_bandwidth=numpy.array([1e6]),
423+
phasecentre=None, kind='cubic', applybeam=True, flux_limit=0.1,
424+
flux_max=numpy.inf, flux_threshold=1.0, insert_method='Nearest',
425+
telescope='LOW') -> SkyModel:
426+
"""Create LOW test skymodel from the GLEAM survey
427+
428+
Stokes I is estimated from a cubic spline fit to the measured fluxes. The polarised flux is always zero.
429+
430+
See http://www.mwatelescope.org/science/gleam-survey The catalog is available from Vizier.
431+
432+
VIII/100 GaLactic and Extragalactic All-sky MWA survey (Hurley-Walker+, 2016)
433+
434+
GaLactic and Extragalactic All-sky Murchison Wide Field Array (GLEAM) survey. I: A low-frequency extragalactic
435+
catalogue. Hurley-Walker N., et al., Mon. Not. R. Astron. Soc., 464, 1146-1167 (2017), 2017MNRAS.464.1146H
436+
437+
:param telescope:
438+
:param npixel: Number of pixels
439+
:param polarisation_frame: Polarisation frame (default PolarisationFrame("stokesI"))
440+
:param cellsize: cellsize in radians
441+
:param frequency:
442+
:param channel_bandwidth: Channel width (Hz)
443+
:param phasecentre: phasecentre (SkyCoord)
444+
:param kind: Kind of interpolation (see scipy.interpolate.interp1d) Default: cubic
445+
:param applybeam: Apply the primary beam?
446+
:param flux_limit: Weakest component
447+
:param flux_max: Maximum strength component to be included in components
448+
:param flux_threshold: Split between components (brighter) and image (weaker)
449+
:param insert_method: Nearest | PSWF | Lanczos
450+
:return:
451+
:return: SkyModel
452+
453+
"""
454+
455+
if phasecentre is None:
456+
phasecentre = SkyCoord(ra=+15.0 * u.deg, dec=-35.0 * u.deg, frame='icrs', equinox='J2000')
457+
458+
radius = npixel * cellsize
459+
460+
sc = create_low_test_skycomponents_from_gleam(flux_limit=flux_limit, polarisation_frame=polarisation_frame,
461+
frequency=frequency, phasecentre=phasecentre,
462+
kind=kind, radius=radius)
463+
464+
sc = filter_skycomponents_by_flux(sc, flux_max=flux_max)
465+
if polarisation_frame is None:
466+
polarisation_frame = PolarisationFrame("stokesI")
467+
468+
npol = polarisation_frame.npol
469+
nchan = len(frequency)
470+
shape = [nchan, npol, npixel, npixel]
471+
w = WCS(naxis=4)
472+
# The negation in the longitude is needed by definition of RA, DEC
473+
w.wcs.cdelt = [-cellsize * 180.0 / numpy.pi, cellsize * 180.0 / numpy.pi, 1.0, channel_bandwidth[0]]
474+
w.wcs.crpix = [npixel // 2 + 1, npixel // 2 + 1, 1.0, 1.0]
475+
w.wcs.ctype = ["RA---SIN", "DEC--SIN", 'STOKES', 'FREQ']
476+
w.wcs.crval = [phasecentre.ra.deg, phasecentre.dec.deg, 1.0, frequency[0]]
477+
w.naxis = 4
478+
w.wcs.radesys = 'ICRS'
479+
w.wcs.equinox = 2000.0
480+
481+
model = create_image_from_array(numpy.zeros(shape), w, polarisation_frame=polarisation_frame)
482+
483+
if applybeam:
484+
beam = create_pb(model, telescope=telescope)
485+
sc = apply_beam_to_skycomponent(sc, beam)
486+
487+
weaksc = filter_skycomponents_by_flux(sc, flux_max=flux_threshold)
488+
brightsc = filter_skycomponents_by_flux(sc, flux_min=flux_threshold)
489+
model = insert_skycomponent(model, weaksc, insert_method=insert_method)
490+
log.info('create_low_test_skymodel_from_gleam: %d bright sources above flux threshold %.3f, %d weak sources below ' %
491+
(len(brightsc), flux_threshold, len(weaksc)))
492+
493+
return SkyModel(components=brightsc, images=[model])
494+
495+
421496
def create_low_test_skycomponents_from_gleam(flux_limit=0.1, polarisation_frame=PolarisationFrame("stokesI"),
422497
frequency=numpy.array([1e8]), kind='cubic', phasecentre=None,
423498
radius=1.0) \
@@ -434,7 +509,6 @@ def create_low_test_skycomponents_from_gleam(flux_limit=0.1, polarisation_frame=
434509
catalogue. Hurley-Walker N., et al., Mon. Not. R. Astron. Soc., 464, 1146-1167 (2017), 2017MNRAS.464.1146H
435510
436511
437-
:rtype: Union[None, List[processing_library.data_models.data_models.Skycomponent], List]
438512
:param flux_limit: Only write components brighter than this (Jy)
439513
:param polarisation_frame: Polarisation frame (default PolarisationFrame("stokesI"))
440514
:param frequency: Frequencies at which the flux will be estimated
@@ -487,13 +561,13 @@ def create_low_test_skycomponents_from_gleam(flux_limit=0.1, polarisation_frame=
487561
gleam_flux_freq[:, i] = filtered_recs['int_flux_%03d' % (f)][:]
488562

489563
skycomps = []
490-
564+
491565
directions = SkyCoord(ra=ras * u.deg, dec=decs * u.deg)
492566
if phasecentre is not None:
493567
separations = directions.separation(phasecentre).to('rad').value
494568
else:
495569
separations = numpy.zeros(len(names))
496-
570+
497571
for isource, name in enumerate(names):
498572
direction = directions[isource]
499573
if separations[isource] < radius:
@@ -692,15 +766,14 @@ def create_unittest_components(model, flux, applypb=False, telescope='LOW', npix
692766
spacing_pixels = int(scale * npixel) // 4
693767
log.info('Spacing in pixels = %s' % spacing_pixels)
694768

695-
696769
if not symmetric:
697770
centers = [(0.2, 1.1)]
698771
else:
699772
centers = list()
700-
773+
701774
if not single:
702775
centers.append([0.0, 0.0])
703-
776+
704777
for x in numpy.linspace(-1.2, 1.2, 7):
705778
if abs(x) > 1e-15:
706779
centers.append([x, x])

0 commit comments

Comments
 (0)