Open in Colab: https://colab.research.google.com/github/casangi/cngi_prototype/blob/master/docs/imaging/mosaic_image_example.ipynb

Mosaic Imaging

This notebook will demonstrate how to create a joint mosaic image. Data is taken from all the pointings in an input visibility dataset and combined to have a single phase-reference center.

This walkthrough is designed to be run in a Jupyter notebook on Google Colaboratory. To open the notebook in colab, go here.

Installation

[1]:

import os
os.system("pip install cngi-prototype==0.0.91")
print('complete')

complete

Dataset

The simulated dataset consists of three fields which contain four point sources over three frequency channels. The ALMA layout alma.cycle6.3.cfg is used, which can be found here.

[2]:

!gdown -q --id 1KzWk0Xg8-xpljTL6m8WEE8KbfLTpskRg
!unzip alma12m_3field_dovpTrue.vis.zarr.zip > /dev/null

!gdown -q --id 1YdJGBi2qtdCuJ6dm4xrXU5zvXeJc7w3v
!unzip alma12m_3field_dovpTrue_gridder_mosaic.img.zarr.zip > /dev/null

#%matplotlib widget

Load Dataset

[3]:

import xarray as xr
from cngi.dio import read_vis

xr.set_options(display_style="html")

infile = "alma12m_3field_dovpTrue.vis.zarr"
mxds = read_vis(infile)
mxds

overwrite_encoded_chunks True

[3]:

<xarray.Dataset>
Dimensions:           (antenna_ids: 43, feed_ids: 43, field_ids: 3, observation_ids: 1, polarization_ids: 1, source_ids: 3, spw_ids: 1, state_ids: 1)
Coordinates:
  * antenna_ids       (antenna_ids) int64 0 1 2 3 4 5 6 ... 36 37 38 39 40 41 42
    antennas          (antenna_ids) <U16 'A001' 'A002' 'A007' ... 'A085' 'A088'
  * field_ids         (field_ids) int64 0 1 2
    fields            (field_ids) <U16 'field_1' 'field_2' 'field_3'
  * feed_ids          (feed_ids) int32 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0
  * observation_ids   (observation_ids) int64 0
    observations      (observation_ids) <U16 'CASA simulation'
  * polarization_ids  (polarization_ids) int64 0
  * source_ids        (source_ids) int32 0 1 2
    sources           (source_ids) <U16 'field_1' 'field_2' 'field_3'
  * spw_ids           (spw_ids) int64 0
  * state_ids         (state_ids) int64 0
Data variables:
    *empty*
Attributes:
    xds0:             <xarray.Dataset>\nDimensions:                (baseline:...
    ANTENNA:          <xarray.Dataset>\nDimensions:        (antenna_id: 43, d...
    FEED:             <xarray.Dataset>\nDimensions:             (d0: 43, d1: ...
    FIELD:            <xarray.Dataset>\nDimensions:        (d1: 1, d2: 2, fie...
    OBSERVATION:      <xarray.Dataset>\nDimensions:         (d1: 2, observati...
    POINTING:         <xarray.Dataset>\nDimensions:      (antenna_id: 43, d2:...
    POLARIZATION:     <xarray.Dataset>\nDimensions:       (d0: 1, d1: 1, d2: ...
    SOURCE:           <xarray.Dataset>\nDimensions:             (d0: 3, d1: 2...
    SPECTRAL_WINDOW:  <xarray.Dataset>\nDimensions:             (d1: 3, spect...
    STATE:            <xarray.Dataset>\nDimensions:   (state_id: 1)\nCoordina...
[4]:

mxds.xds0

[4]:

<xarray.Dataset>
Dimensions:                (baseline: 903, chan: 3, pol: 1, pol_id: 1, spw_id: 1, time: 192, uvw_index: 3)
Coordinates:
  * baseline               (baseline) int64 0 1 2 3 4 5 ... 898 899 900 901 902
  * chan                   (chan) float64 3.4e+11 3.74e+11 4.08e+11
    chan_width             (chan) float64 dask.array<chunksize=(1,), meta=np.ndarray>
    effective_bw           (chan) float64 dask.array<chunksize=(1,), meta=np.ndarray>
  * pol                    (pol) int32 9
  * pol_id                 (pol_id) int32 0
    resolution             (chan) float64 dask.array<chunksize=(1,), meta=np.ndarray>
  * spw_id                 (spw_id) int32 0
  * time                   (time) datetime64[ns] 2011-05-27T04:32:14.77260398...
Dimensions without coordinates: uvw_index
Data variables: (12/20)
    ANTENNA1               (baseline) int32 dask.array<chunksize=(301,), meta=np.ndarray>
    ANTENNA2               (baseline) int32 dask.array<chunksize=(301,), meta=np.ndarray>
    ARRAY_ID               (time, baseline) int32 dask.array<chunksize=(32, 301), meta=np.ndarray>
    CORRECTED_DATA         (time, baseline, chan, pol) complex128 dask.array<chunksize=(32, 301, 1, 1), meta=np.ndarray>
    CORRECTED_DATA_WEIGHT  (time, baseline, chan, pol) float64 dask.array<chunksize=(32, 301, 1, 1), meta=np.ndarray>
    DATA                   (time, baseline, chan, pol) complex128 dask.array<chunksize=(32, 301, 1, 1), meta=np.ndarray>
    ...                     ...
    OBSERVATION_ID         (time, baseline) int32 dask.array<chunksize=(32, 301), meta=np.ndarray>
    PROCESSOR_ID           (time, baseline) int32 dask.array<chunksize=(32, 301), meta=np.ndarray>
    SCAN_NUMBER            (time, baseline) int32 dask.array<chunksize=(32, 301), meta=np.ndarray>
    STATE_ID               (time, baseline) int32 dask.array<chunksize=(32, 301), meta=np.ndarray>
    TIME_CENTROID          (time, baseline) float64 dask.array<chunksize=(32, 301), meta=np.ndarray>
    UVW                    (time, baseline, uvw_index) float64 dask.array<chunksize=(32, 301, 3), meta=np.ndarray>
Attributes:
    corr_product:     [[0, 0]]
    data_groups:      [{'0': {'data': 'DATA', 'flag': 'FLAG', 'id': '0', 'uvw...
    freq_group:       0
    freq_group_name:  Group 1
    if_conv_chain:    0
    meas_freq_ref:    1
    name:             Band_8
    net_sideband:     1
    num_chan:         3
    num_corr:         1
    ref_frequency:    340000000000.0
    total_bandwidth:  101999999999.99998

Grid Parameters

[5]:

grid_parms = {}
grid_parms['chan_mode'] = 'cube'
grid_parms['image_size'] = [1000,720]
grid_parms['cell_size'] = [0.04,0.04]
grid_parms['fft_padding'] = 1.0
grid_parms['phase_center'] = mxds.FIELD.PHASE_DIR[1,0,:].data.compute()

Direction Rotation

The UVW coordinates must be rotated and the visbility DATA must be phase rotated, relative to the mosaic phase center specified by rotation_parms['image_phase_center'].

direction_rotate documentation

[6]:

from ngcasa.imaging import direction_rotate
import numpy as np
import dask

xr.set_options(display_style="html")

infile = "alma12m_3field_dovpTrue.vis.zarr"
#DATA shape (192, 903, 3, 1), ZARR chunking (32, 301, 1, 1)
mxds = read_vis(infile,chunks={'time':192,'baseline':903,'chan':1})

sel_parms = {}
sel_parms['xds'] = 'xds0'

rotation_parms = {}
rotation_parms['new_phase_center'] = grid_parms['phase_center']
rotation_parms['common_tangent_reprojection'] = True
rotation_parms['single_precision'] = False

mxds = direction_rotate(mxds, rotation_parms, sel_parms)

mxds.xds0

overwrite_encoded_chunks True
######################### Start direction_rotate #########################
Setting default data_group_in  to  {'data': 'DATA', 'flag': 'FLAG', 'id': '0', 'uvw': 'UVW', 'weight': 'DATA_WEIGHT'}
Setting default data_group_out  to  {'data': 'DATA_ROT', 'flag': 'FLAG', 'id': '2', 'uvw': 'UVW_ROT', 'weight': 'DATA_WEIGHT'}
######################### Created graph for direction_rotate #########################

[6]:

<xarray.Dataset>
Dimensions:                (baseline: 903, chan: 3, pol: 1, pol_id: 1, spw_id: 1, time: 192, uvw_index: 3)
Coordinates:
  * baseline               (baseline) int64 0 1 2 3 4 5 ... 898 899 900 901 902
  * chan                   (chan) float64 3.4e+11 3.74e+11 4.08e+11
  * pol                    (pol) int32 9
  * time                   (time) datetime64[ns] 2011-05-27T04:32:14.77260398...
    chan_width             (chan) float64 dask.array<chunksize=(1,), meta=np.ndarray>
    effective_bw           (chan) float64 dask.array<chunksize=(1,), meta=np.ndarray>
  * pol_id                 (pol_id) int32 0
    resolution             (chan) float64 dask.array<chunksize=(1,), meta=np.ndarray>
  * spw_id                 (spw_id) int32 0
Dimensions without coordinates: uvw_index
Data variables: (12/22)
    ANTENNA1               (baseline) int32 dask.array<chunksize=(903,), meta=np.ndarray>
    ANTENNA2               (baseline) int32 dask.array<chunksize=(903,), meta=np.ndarray>
    ARRAY_ID               (time, baseline) int32 dask.array<chunksize=(192, 903), meta=np.ndarray>
    CORRECTED_DATA         (time, baseline, chan, pol) complex128 dask.array<chunksize=(192, 903, 1, 1), meta=np.ndarray>
    CORRECTED_DATA_WEIGHT  (time, baseline, chan, pol) float64 dask.array<chunksize=(192, 903, 1, 1), meta=np.ndarray>
    DATA                   (time, baseline, chan, pol) complex128 dask.array<chunksize=(192, 903, 1, 1), meta=np.ndarray>
    ...                     ...
    SCAN_NUMBER            (time, baseline) int32 dask.array<chunksize=(192, 903), meta=np.ndarray>
    STATE_ID               (time, baseline) int32 dask.array<chunksize=(192, 903), meta=np.ndarray>
    TIME_CENTROID          (time, baseline) float64 dask.array<chunksize=(192, 903), meta=np.ndarray>
    UVW                    (time, baseline, uvw_index) float64 dask.array<chunksize=(192, 903, 3), meta=np.ndarray>
    UVW_ROT                (time, baseline, uvw_index) float64 dask.array<chunksize=(192, 903, 3), meta=np.ndarray>
    DATA_ROT               (time, baseline, chan, pol) complex128 dask.array<chunksize=(192, 903, 1, 1), meta=np.ndarray>
Attributes:
    corr_product:     [[0, 0]]
    data_groups:      [{'0': {'data': 'DATA', 'flag': 'FLAG', 'id': '0', 'uvw...
    freq_group:       0
    freq_group_name:  Group 1
    if_conv_chain:    0
    meas_freq_ref:    1
    name:             Band_8
    net_sideband:     1
    num_chan:         3
    num_corr:         1
    ref_frequency:    340000000000.0
    total_bandwidth:  101999999999.99998

Make Imaging Weights

make_imaging_weight documentation

[7]:

from ngcasa.imaging import make_imaging_weight

imaging_weights_parms = {}
imaging_weights_parms['weighting'] = 'natural'

sel_parms = {}
sel_parms['xds'] = 'xds0'
sel_parms['data_group_in_id'] = 2

mxds = make_imaging_weight(mxds, imaging_weights_parms, grid_parms, sel_parms)

imaging_weights_parms = {}
imaging_weights_parms['weighting'] = 'natural'

sel_parms = {}
sel_parms['xds'] = 'xds0'
sel_parms['data_group_in_id'] = 0

mxds = make_imaging_weight(mxds, imaging_weights_parms, grid_parms, sel_parms)
mxds.xds0

######################### Start make_imaging_weights #########################
Setting data_group_in  to  {'data': 'DATA_ROT', 'flag': 'FLAG', 'id': '2', 'uvw': 'UVW_ROT', 'weight': 'DATA_WEIGHT'}
Setting default data_group_out  to  {'data': 'DATA_ROT', 'flag': 'FLAG', 'id': '2', 'uvw': 'UVW_ROT', 'weight': 'DATA_WEIGHT', 'imaging_weight': 'IMAGING_WEIGHT'}
Since weighting is natural input weight will be reused as imaging weight.
######################### Created graph for make_imaging_weight #########################
######################### Start make_imaging_weights #########################
Setting data_group_in  to  {'data': 'DATA', 'flag': 'FLAG', 'id': '0', 'uvw': 'UVW', 'weight': 'DATA_WEIGHT'}
Setting default data_group_out  to  {'data': 'DATA', 'flag': 'FLAG', 'id': '0', 'uvw': 'UVW', 'weight': 'DATA_WEIGHT', 'imaging_weight': 'IMAGING_WEIGHT'}
Since weighting is natural input weight will be reused as imaging weight.
######################### Created graph for make_imaging_weight #########################

[7]:

<xarray.Dataset>
Dimensions:                (baseline: 903, chan: 3, pol: 1, pol_id: 1, spw_id: 1, time: 192, uvw_index: 3)
Coordinates:
  * baseline               (baseline) int64 0 1 2 3 4 5 ... 898 899 900 901 902
  * chan                   (chan) float64 3.4e+11 3.74e+11 4.08e+11
  * pol                    (pol) int32 9
  * time                   (time) datetime64[ns] 2011-05-27T04:32:14.77260398...
    chan_width             (chan) float64 dask.array<chunksize=(1,), meta=np.ndarray>
    effective_bw           (chan) float64 dask.array<chunksize=(1,), meta=np.ndarray>
  * pol_id                 (pol_id) int32 0
    resolution             (chan) float64 dask.array<chunksize=(1,), meta=np.ndarray>
  * spw_id                 (spw_id) int32 0
Dimensions without coordinates: uvw_index
Data variables: (12/22)
    ANTENNA1               (baseline) int32 dask.array<chunksize=(903,), meta=np.ndarray>
    ANTENNA2               (baseline) int32 dask.array<chunksize=(903,), meta=np.ndarray>
    ARRAY_ID               (time, baseline) int32 dask.array<chunksize=(192, 903), meta=np.ndarray>
    CORRECTED_DATA         (time, baseline, chan, pol) complex128 dask.array<chunksize=(192, 903, 1, 1), meta=np.ndarray>
    CORRECTED_DATA_WEIGHT  (time, baseline, chan, pol) float64 dask.array<chunksize=(192, 903, 1, 1), meta=np.ndarray>
    DATA                   (time, baseline, chan, pol) complex128 dask.array<chunksize=(192, 903, 1, 1), meta=np.ndarray>
    ...                     ...
    SCAN_NUMBER            (time, baseline) int32 dask.array<chunksize=(192, 903), meta=np.ndarray>
    STATE_ID               (time, baseline) int32 dask.array<chunksize=(192, 903), meta=np.ndarray>
    TIME_CENTROID          (time, baseline) float64 dask.array<chunksize=(192, 903), meta=np.ndarray>
    UVW                    (time, baseline, uvw_index) float64 dask.array<chunksize=(192, 903, 3), meta=np.ndarray>
    UVW_ROT                (time, baseline, uvw_index) float64 dask.array<chunksize=(192, 903, 3), meta=np.ndarray>
    DATA_ROT               (time, baseline, chan, pol) complex128 dask.array<chunksize=(192, 903, 1, 1), meta=np.ndarray>
Attributes:
    corr_product:     [[0, 0]]
    data_groups:      [{'0': {'data': 'DATA', 'flag': 'FLAG', 'id': '0', 'uvw...
    freq_group:       0
    freq_group_name:  Group 1
    if_conv_chain:    0
    meas_freq_ref:    1
    name:             Band_8
    net_sideband:     1
    num_chan:         3
    num_corr:         1
    ref_frequency:    340000000000.0
    total_bandwidth:  101999999999.99998

Make Gridding Convolution Functions

make_gridding_convolution_function

[8]:

from ngcasa.imaging import make_gridding_convolution_function
import numpy as np
import dask.array as da
from cngi.dio import write_image

gcf_parms = {}
gcf_parms['function'] = 'casa_airy'
gcf_parms['list_dish_diameters'] = np.array([10.7])
gcf_parms['list_blockage_diameters'] = np.array([0.75])

unique_ant_indx = mxds.ANTENNA.DISH_DIAMETER.values
unique_ant_indx[unique_ant_indx == 12.0] = 0

gcf_parms['unique_ant_indx'] = unique_ant_indx.astype(int)
gcf_parms['phase_center'] = grid_parms['phase_center']

sel_parms = {}
sel_parms['xds'] = 'xds0'
sel_parms['data_group_in_id'] = 2

gcf_xds = make_gridding_convolution_function(mxds, gcf_parms, grid_parms, sel_parms)
write_image(gcf_xds,'mosaic_gcf.gcf.zarr')
gcf_xds = xr.open_zarr('mosaic_gcf.gcf.zarr')
gcf_xds

######################### Start make_gridding_convolution_function #########################
Setting data_group_in  to  {'data': 'DATA_ROT', 'flag': 'FLAG', 'id': '2', 'uvw': 'UVW_ROT', 'weight': 'DATA_WEIGHT', 'imaging_weight': 'DATA_WEIGHT'}
Setting default data_group_out  to  {'data': 'DATA_ROT', 'flag': 'FLAG', 'id': '3', 'uvw': 'UVW_ROT', 'weight': 'DATA_WEIGHT', 'imaging_weight': 'DATA_WEIGHT'}
Setting default chan_tolerance_factor  to  0.005
Setting default oversampling  to  [10, 10]
Setting default max_support  to  [15, 15]
Setting default support_cut_level  to  0.025
Setting default a_chan_num_chunk  to  3
Setting default image_center  to  [500 360]
#########################  Created graph for make_gridding_convolution_function #########################
Time to store and execute graph  write_zarr 7.4372382164001465

[8]:

<xarray.Dataset>
Dimensions:             (baseline: 903, chan: 3, conv_baseline: 1, conv_chan: 3, conv_pol: 1, field_id: 3, l: 1000, m: 720, pol: 1, u: 160, v: 160, xy: 2)
Coordinates:
  * field_id            (field_id) int64 0 1 2
  * l                   (l) int64 0 1 2 3 4 5 6 ... 993 994 995 996 997 998 999
  * m                   (m) int64 0 1 2 3 4 5 6 ... 713 714 715 716 717 718 719
  * u                   (u) int64 0 1 2 3 4 5 6 ... 153 154 155 156 157 158 159
  * v                   (v) int64 0 1 2 3 4 5 6 ... 153 154 155 156 157 158 159
  * xy                  (xy) int64 0 1
Dimensions without coordinates: baseline, chan, conv_baseline, conv_chan, conv_pol, pol
Data variables:
    CF_BASELINE_MAP     (baseline) int64 dask.array<chunksize=(903,), meta=np.ndarray>
    CF_CHAN_MAP         (chan) int64 dask.array<chunksize=(1,), meta=np.ndarray>
    CF_POL_MAP          (pol) int64 dask.array<chunksize=(1,), meta=np.ndarray>
    CONV_KERNEL         (conv_baseline, conv_chan, conv_pol, u, v) float64 dask.array<chunksize=(1, 1, 1, 160, 160), meta=np.ndarray>
    PHASE_GRADIENT      (field_id, u, v) complex128 dask.array<chunksize=(1, 160, 160), meta=np.ndarray>
    PS_CORR_IMAGE       (l, m) float64 dask.array<chunksize=(1000, 720), meta=np.ndarray>
    SUPPORT             (conv_baseline, conv_chan, conv_pol, xy) int64 dask.array<chunksize=(1, 1, 1, 2), meta=np.ndarray>
    WEIGHT_CONV_KERNEL  (conv_baseline, conv_chan, conv_pol, u, v) float64 dask.array<chunksize=(1, 1, 1, 160, 160), meta=np.ndarray>
Attributes:
    cell_uv:          [-515.6620156177408, 716.197243913529]
    oversampling:     [10, 10]
    write_zarr_time:  7.4372382164001465

Make Mosaic Primary Beam, PSF, and Image

make_mosaic_pb

make_psf

make_image_with_gcf

[9]:

from ngcasa.imaging import make_mosaic_pb
from cngi.dio import read_image

vis_sel_parms = {}
vis_sel_parms['xds'] = 'xds0'
vis_sel_parms['data_group_in_id'] = 2

img_sel_parms = {}
img_xds= xr.Dataset()

img_xds = make_mosaic_pb(mxds,gcf_xds,img_xds,vis_sel_parms,img_sel_parms,grid_parms)

###############################################

from ngcasa.imaging import make_psf

vis_sel_parms = {}
vis_sel_parms['xds'] = 'xds0'
vis_sel_parms['data_group_in_id'] = 2

img_sel_parms = {}
img_sel_parms['data_group_out_id'] = 0

img_xds = make_psf(mxds, img_xds, grid_parms, vis_sel_parms, img_sel_parms)

##############################################

from ngcasa.imaging import make_image_with_gcf
from cngi.dio import write_image

vis_select_parms = {}
vis_select_parms['xds'] = 'xds0'
vis_select_parms['data_group_in_id'] = 2

img_select_parms = {}
img_select_parms['data_group_in_id'] = 0
img_select_parms['data_group_out_id'] = 0

norm_parms = {}
norm_parms['norm_type'] = 'flat_sky'

img_xds = make_image_with_gcf(mxds,gcf_xds, img_xds, grid_parms, norm_parms, vis_select_parms, img_select_parms)
write_image(img_xds,'mosaic_img.img.zarr')
img_xds = read_image('mosaic_img.img.zarr')
img_xds

######################### Start make_mosaic_pb #########################
Setting default image_center  to  [500 360]
Setting data_group_in  to  {'data': 'DATA_ROT', 'flag': 'FLAG', 'id': '2', 'uvw': 'UVW_ROT', 'weight': 'DATA_WEIGHT', 'imaging_weight': 'DATA_WEIGHT'}
Setting default data_group_out  to  {'data': 'DATA_ROT', 'flag': 'FLAG', 'id': '3', 'uvw': 'UVW_ROT', 'weight': 'DATA_WEIGHT', 'imaging_weight': 'DATA_WEIGHT'}
Setting default data_group_in  to  {'id': '0'}
Setting default data_group_out  to  {'id': '0', 'pb': 'PB', 'weight_pb': 'WEIGHT_PB', 'weight_pb_sum_weight': 'WEIGHT_PB_SUM_WEIGHT'}
#########################  Created graph for make_mosaic_pb #########################
######################### Start make_psf #########################
Setting default image_center  to  [500 360]
Setting data_group_in  to  {'data': 'DATA_ROT', 'flag': 'FLAG', 'id': '2', 'uvw': 'UVW_ROT', 'weight': 'DATA_WEIGHT', 'imaging_weight': 'DATA_WEIGHT'}
Setting default data_group_out  to  {'data': 'DATA_ROT', 'flag': 'FLAG', 'id': '3', 'uvw': 'UVW_ROT', 'weight': 'DATA_WEIGHT', 'imaging_weight': 'DATA_WEIGHT'}
Setting default data_group_in  to  {'id': '0', 'pb': 'PB', 'weight_pb': 'WEIGHT_PB', 'weight_pb_sum_weight': 'WEIGHT_PB_SUM_WEIGHT'}
Setting default data_group_out [' pb ']  to  PB
Setting default data_group_out [' weight_pb ']  to  WEIGHT_PB
Setting default data_group_out [' weight_pb_sum_weight ']  to  WEIGHT_PB_SUM_WEIGHT
Setting default data_group_out [' sum_weight ']  to  PSF_SUM_WEIGHT
Setting default data_group_out [' psf ']  to  PSF
Setting default data_group_out [' psf_fit ']  to  PSF_FIT
######################### Created graph for make_psf #########################
######################### Start make_image_with_gcf #########################
Setting default image_center  to  [500 360]
Setting default single_precision  to  True
Setting default pb_limit  to  0.2
Setting data_group_in  to  {'data': 'DATA_ROT', 'flag': 'FLAG', 'id': '2', 'uvw': 'UVW_ROT', 'weight': 'DATA_WEIGHT', 'imaging_weight': 'DATA_WEIGHT'}
Setting default data_group_out  to  {'data': 'DATA_ROT', 'flag': 'FLAG', 'id': '3', 'uvw': 'UVW_ROT', 'weight': 'DATA_WEIGHT', 'imaging_weight': 'DATA_WEIGHT'}
Setting data_group_in  to  {'id': '0', 'pb': 'PB', 'weight_pb': 'WEIGHT_PB', 'weight_pb_sum_weight': 'WEIGHT_PB_SUM_WEIGHT', 'sum_weight': 'PSF_SUM_WEIGHT', 'psf': 'PSF', 'psf_fit': 'PSF_FIT'}
Setting default data_group_out [' pb ']  to  PB
Setting default data_group_out [' weight_pb ']  to  WEIGHT_PB
Setting default data_group_out [' weight_pb_sum_weight ']  to  WEIGHT_PB_SUM_WEIGHT
Setting default data_group_out [' sum_weight ']  to  SUM_WEIGHT
Setting default data_group_out [' psf ']  to  PSF
Setting default data_group_out [' psf_fit ']  to  PSF_FIT
Setting default data_group_out [' image ']  to  IMAGE
grid sizes 1000 720
#########################  Created graph for make_mosaic_with_gcf #########################
Time to store and execute graph  write_zarr 10.509612560272217

[9]:

<xarray.Dataset>
Dimensions:               (chan: 3, elps_index: 3, l: 1000, m: 720, pol: 1, time: 1)
Coordinates:
  * chan                  (chan) float64 3.4e+11 3.74e+11 4.08e+11
    chan_width            (chan) float64 dask.array<chunksize=(1,), meta=np.ndarray>
    declination           (l, m) float64 dask.array<chunksize=(250, 180), meta=np.ndarray>
  * l                     (l) float64 9.696e-05 9.677e-05 ... -9.677e-05
  * m                     (m) float64 -6.981e-05 -6.962e-05 ... 6.962e-05
  * pol                   (pol) int32 9
    right_ascension       (l, m) float64 dask.array<chunksize=(250, 180), meta=np.ndarray>
  * time                  (time) datetime64[ns] 2011-05-27T20:27:14.772603989
Dimensions without coordinates: elps_index
Data variables:
    IMAGE                 (l, m, time, chan, pol) float64 dask.array<chunksize=(1000, 720, 1, 1, 1), meta=np.ndarray>
    PB                    (l, m, time, chan, pol) float64 dask.array<chunksize=(1000, 720, 1, 1, 1), meta=np.ndarray>
    PSF                   (l, m, time, chan, pol) float64 dask.array<chunksize=(1000, 720, 1, 1, 1), meta=np.ndarray>
    PSF_FIT               (time, chan, pol, elps_index) float64 dask.array<chunksize=(1, 1, 1, 3), meta=np.ndarray>
    PSF_SUM_WEIGHT        (time, chan, pol) float64 dask.array<chunksize=(1, 1, 1), meta=np.ndarray>
    SUM_WEIGHT            (time, chan, pol) float64 dask.array<chunksize=(1, 1, 1), meta=np.ndarray>
    WEIGHT_PB             (l, m, time, chan, pol) float64 dask.array<chunksize=(1000, 720, 1, 1, 1), meta=np.ndarray>
    WEIGHT_PB_SUM_WEIGHT  (time, chan, pol) float64 dask.array<chunksize=(1, 1, 1), meta=np.ndarray>
Attributes:
    axis_units:           ['rad', 'rad', 'time', 'Hz', 'pol']
    data_groups:          [{'0': {'id': '0', 'image': 'IMAGE', 'pb': 'PB', 'p...
    direction_reference:  FK5
    spectral_reference:   lsrk
    unit:                 Jy/beam
    velocity_type:        radio
    write_zarr_time:      10.509612560272217

Compare CASA and ngCASA Primary Beams

[10]:

import matplotlib.pylab as plt
import numpy as np
from ipywidgets import interactive
import scipy
from scipy.signal import decimate
from cngi.image import implot

img_xds = read_image('mosaic_img.img.zarr').isel(time=0,pol=0)
casa_img_xds = read_image('alma12m_3field_dovpTrue_gridder_mosaic.img.zarr').isel(time=0,pol=0)
pb_limit = 0.2
extent = extent=(np.min(casa_img_xds.m),np.max(casa_img_xds.m),np.min(casa_img_xds.l),np.max(casa_img_xds.l))

def comparison_plots(chan):
    plt.close('all')
    print('Frequency',img_xds.chan[chan].values/10**9, 'GHz')
    mosaic_pb = img_xds.PB.isel(chan=chan)
    mosaic_pb = mosaic_pb.where(mosaic_pb > pb_limit,other=np.nan)

    casa_mosaic_pb = casa_img_xds.PB.isel(chan=chan)
    casa_mosaic_pb = casa_mosaic_pb.where(casa_mosaic_pb > pb_limit,other=np.nan)

    fig0, ax0 = plt.subplots(1, 2, sharey=True,figsize=(10, 5))
    im0 = ax0[0].imshow(mosaic_pb,extent=extent,cmap='jet')
    im1 = ax0[1].imshow(casa_mosaic_pb,extent=extent,cmap='jet')
    ax0[0].title.set_text('ngCASA Mosaic PB')
    ax0[1].title.set_text('CASA Mosaic PB')
    ax0[0].set_xlabel('m'), ax0[1].set_xlabel('m'), ax0[0].set_ylabel('l'), ax0[1].set_ylabel('l')
    fig0.colorbar(im0, ax=ax0[0], fraction=0.046, pad=0.04)
    fig0.colorbar(im1, ax=ax0[1], fraction=0.046, pad=0.04)

    plt.figure()
    plt.plot(casa_mosaic_pb.l,mosaic_pb.isel(m=360),label='ngCASA PB')
    plt.plot(casa_mosaic_pb.l,casa_mosaic_pb.isel(m=360),'*',label='CASA PB',markersize=1)
    plt.legend()
    plt.xlabel('l')
    plt.ylabel('Amplitude')
    plt.title('Mosaic PB Cross Section')

    diff_image = mosaic_pb - casa_mosaic_pb

    plt.figure()
    plt.imshow(100*(mosaic_pb - casa_mosaic_pb),extent=extent,cmap='jet')
    plt.xlabel('m'), plt.ylabel('l')
    plt.colorbar()
    plt.title('Percentage Difference Mosaic PB')

    plt.show()

interactive_plot = interactive(comparison_plots, chan=(0, 2))
output = interactive_plot.children[-1]
output.layout.auto_scroll_threshold = 9999;
#interactive_plot
comparison_plots(1)


Frequency 374.0 GHz
../_images/imaging_mosaic_image_example_19_1.png
../_images/imaging_mosaic_image_example_19_2.png
../_images/imaging_mosaic_image_example_19_3.png

Get Simulated Sources l,m Coordinates

[11]:

from astropy.coordinates import SkyCoord
from astropy.wcs import WCS
rad_to_deg =  180/np.pi
deg_to_rad = np.pi/180
arcsec_to_deg = 1/3600
arcsec_to_rad = np.pi/(180*3600)

phase_center = grid_parms['phase_center']
w = WCS(naxis=2)
w.wcs.crpix = np.array(grid_parms['image_size'])//2
w.wcs.cdelt = np.array([-0.04,0.04])*arcsec_to_deg
w.wcs.crval = phase_center*rad_to_deg
w.wcs.ctype = ['RA---SIN','DEC--SIN']

ra = ['12h01m51.903005s','12h01m52.430856s','12h01m52.958707s','12h01m52.259s','12h01m52s','12h01m53.153s']
dec = ['-18d51m49.94373s','-18d51m49.94369s','-18d51m49.94365s','-18d51m42.983s','-18d51m46s','-18d51m59.305s']
ps_skycoord = SkyCoord(ra=ra,dec=dec,frame='fk5')


ra_dec = np.array([ps_skycoord.ra.degree,ps_skycoord.dec.degree]).T
lm_pix_pos = w.all_world2pix(ra_dec, 1)

cell_size = np.array(grid_parms['cell_size'])*arcsec_to_rad
cell_size[0] = -cell_size[0]
image_center = np.array(grid_parms['image_size'])//2
source_lm_pos = lm_pix_pos*cell_size - image_center*cell_size

Compare CASA and ngCASA Sky Images

[12]:

import matplotlib.pylab as plt
import numpy as np
from ipywidgets import interactive
import scipy
from scipy.signal import decimate

img_xds = read_image('mosaic_img.img.zarr',chunks={'l':grid_parms['image_size'][0],'m':grid_parms['image_size'][1]}).isel(time=0,pol=0)
casa_img_xds = read_image('alma12m_3field_dovpTrue_gridder_mosaic.img.zarr',chunks={'l':grid_parms['image_size'][0],'m':grid_parms['image_size'][1]}).isel(time=0,pol=0)
pb_limit = 0.2
extent = extent=(np.min(casa_img_xds.m),np.max(casa_img_xds.m),np.min(casa_img_xds.l),np.max(casa_img_xds.l))

ngcasa_image_name = 'IMAGE'
casa_image_name = 'IMAGE_PBCOR'

def comparison_plots(chan):
    print('Frequency',img_xds.chan[chan].values/10**9, 'GHz')
    mosaic_pb = img_xds.PB.isel(chan=chan)
    casa_mosaic_pb = casa_img_xds.PB.isel(chan=chan)

    mosaic_img = img_xds[ngcasa_image_name].isel(chan=chan)
    mosaic_img = mosaic_img.where(mosaic_pb > pb_limit,other=np.nan)

    casa_mosaic_img = casa_img_xds[casa_image_name].isel(chan=chan)
    casa_mosaic_img = casa_mosaic_img.where(casa_mosaic_pb > pb_limit,other=np.nan)

    sim_sources = np.array([1.5,1.76,2.0,2.0,1.456,1.888])

    print('################## Flux of Point Sources ##################')
    print('Sim    ','ngCASA  ', 'CASA')
    for i,s in enumerate(sim_sources):
        ngcasa_recovered_val = img_xds[ngcasa_image_name].isel(chan=chan).interp(l=source_lm_pos[i,0],m=source_lm_pos[i,1]).values
        casa_recovered_val = casa_mosaic_img.interp(l=source_lm_pos[i,0],m=source_lm_pos[i,1]).values
        print('{0:.3f}  '.format(s),'{0:.4f}  '.format(ngcasa_recovered_val),'{0:.4f}'.format(casa_recovered_val))


    print('############ Percentage Difference Flux to Sim ############')
    print('ngCASA  ', 'CASA')
    for i,s in enumerate(sim_sources):
        ngcasa_recovered_val = img_xds[ngcasa_image_name].isel(chan=chan).interp(l=source_lm_pos[i,0],m=source_lm_pos[i,1]).values
        casa_recovered_val = casa_mosaic_img.interp(l=source_lm_pos[i,0],m=source_lm_pos[i,1]).values
        print('{0:.4f}  '.format(100*(s-ngcasa_recovered_val)/s),'{0:.4f}'.format(100*(s-casa_recovered_val)/s))

    fig0, ax0 = plt.subplots(1, 2, sharey=True,figsize=(10, 5))
    im0 = ax0[0].imshow(mosaic_img,cmap='jet',extent=extent)
    im1 = ax0[1].imshow(casa_mosaic_img,cmap='jet',extent=extent)
    ax0[0].set_xlabel('m'), ax0[1].set_xlabel('m'), ax0[0].set_ylabel('l'), ax0[1].set_ylabel('l')
    ax0[0].title.set_text('ngCASA Mosaic Image')
    ax0[1].title.set_text('CASA Mosaic Image')
    fig0.colorbar(im0, ax=ax0[0], fraction=0.046, pad=0.04)
    fig0.colorbar(im1, ax=ax0[1], fraction=0.046, pad=0.04)

    plt.figure()
    plt.imshow((100*(mosaic_img - casa_mosaic_img)/2),cmap='jet',extent=extent)
    plt.colorbar()
    plt.xlabel('m'), plt.ylabel('l')
    plt.title('Percentage Difference Mosaic Image')

    plt.show()

interactive_plot = interactive(comparison_plots, chan=(0, 2))
output = interactive_plot.children[-1]
output.layout.height = '1450px'
#interactive_plot
comparison_plots(1)


Frequency 374.0 GHz
################## Flux of Point Sources ##################
Sim     ngCASA   CASA
1.500   1.4935   1.4909
1.760   1.7518   1.7458
2.000   1.9950   1.9913
2.000   2.0391   2.0546
1.456   1.4630   1.4686
1.888   1.8954   1.8425
############ Percentage Difference Flux to Sim ############
ngCASA   CASA
0.4325   0.6079
0.4672   0.8043
0.2492   0.4330
-1.9533   -2.7304
-0.4825   -0.8683
-0.3894   2.4083
../_images/imaging_mosaic_image_example_23_1.png
../_images/imaging_mosaic_image_example_23_2.png