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


Alma Mosaic Imaging Test

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

[2]:
!gdown -q --id 16E-YPVLlgfth-2YVOdL8Uk3h3QT-73ck
!unzip Antennae_North_cal_lsrk.vis.zarr.zip > /dev/null

!gdown -q --id 1C_0ZarxgMo8mEdyDqs7IXf4IPlB7yzs9
!unzip Antennae_North_cal_lsrk_cube.img.zarr.zip > /dev/null

#%matplotlib widget

Load Dataset

[3]:
import xarray as xr
from cngi.dio import read_vis
#%matplotlib widget

xr.set_options(display_style="html")

infile = "Antennae_North_cal_lsrk.vis.zarr"
mxds = read_vis(infile)
print(mxds.xds0)

overwrite_encoded_chunks True
<xarray.Dataset>
Dimensions:         (baseline: 86, chan: 166, pol: 1, pol_id: 1, spw_id: 1, time: 1355, uvw_index: 3)
Coordinates:
  * baseline        (baseline) int64 0 1 2 3 4 5 6 7 ... 78 79 80 81 82 83 84 85
  * chan            (chan) float64 3.449e+11 3.449e+11 ... 3.43e+11 3.43e+11
    chan_width      (chan) float64 dask.array<chunksize=(10,), meta=np.ndarray>
    effective_bw    (chan) float64 dask.array<chunksize=(10,), meta=np.ndarray>
  * pol             (pol) int32 9
  * pol_id          (pol_id) int32 0
    resolution      (chan) float64 dask.array<chunksize=(10,), meta=np.ndarray>
  * spw_id          (spw_id) int32 0
  * time            (time) datetime64[ns] 2011-05-28T01:41:17.375999451 ... 2...
Dimensions without coordinates: uvw_index
Data variables: (12/17)
    ANTENNA1        (baseline) int32 dask.array<chunksize=(86,), meta=np.ndarray>
    ANTENNA2        (baseline) int32 dask.array<chunksize=(86,), meta=np.ndarray>
    ARRAY_ID        (time, baseline) int32 dask.array<chunksize=(1355, 86), meta=np.ndarray>
    DATA            (time, baseline, chan, pol) complex128 dask.array<chunksize=(1355, 86, 10, 1), meta=np.ndarray>
    DATA_WEIGHT     (time, baseline, chan, pol) float64 dask.array<chunksize=(1355, 86, 10, 1), meta=np.ndarray>
    EXPOSURE        (time, baseline) float64 dask.array<chunksize=(1355, 86), meta=np.ndarray>
    ...              ...
    OBSERVATION_ID  (time, baseline) int32 dask.array<chunksize=(1355, 86), meta=np.ndarray>
    PROCESSOR_ID    (time, baseline) int32 dask.array<chunksize=(1355, 86), meta=np.ndarray>
    SCAN_NUMBER     (time, baseline) int32 dask.array<chunksize=(1355, 86), meta=np.ndarray>
    STATE_ID        (time, baseline) int32 dask.array<chunksize=(1355, 86), meta=np.ndarray>
    TIME_CENTROID   (time, baseline) float64 dask.array<chunksize=(1355, 86), meta=np.ndarray>
    UVW             (time, baseline, uvw_index) float64 dask.array<chunksize=(1355, 86, 3), meta=np.ndarray>
Attributes: (12/13)
    bbc_no:           1
    corr_product:     [[0, 0]]
    data_groups:      [{'0': {'data': 'DATA', 'flag': 'FLAG', 'id': '0', 'uvw...
    freq_group:       0
    freq_group_name:
    if_conv_chain:    0
    ...               ...
    name:
    net_sideband:     1
    num_chan:         166
    num_corr:         1
    ref_frequency:    344871541760.4001
    total_bandwidth:  1864427170.9199219

Grid Parameters

[4]:
grid_parms = {}
grid_parms['chan_mode'] = 'cube'
grid_parms['image_size'] = [500,500]
grid_parms['cell_size'] = [0.13,0.13]
grid_parms['fft_padding'] = 1.0
grid_parms['phase_center'] = mxds.FIELD.PHASE_DIR[12,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

[5]:
#Reload modules
%load_ext autoreload
%autoreload 2
[6]:
from ngcasa.imaging import direction_rotate
import numpy as np
import dask
from cngi.vis import apply_flags

mxds0 = apply_flags(mxds, 'xds0', flags='FLAG')

xr.set_options(display_style="html")

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

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


mxds1 = direction_rotate(mxds0, rotation_parms, sel_parms)
mxds1.xds0.data_groups
######################### Start direction_rotate #########################
Setting 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': '1', 'uvw': 'UVW_ROT', 'weight': 'DATA_WEIGHT'}
######################### Created graph for direction_rotate #########################
[6]:
[{'0': {'data': 'DATA',
   'flag': 'FLAG',
   'id': '0',
   'uvw': 'UVW',
   'weight': 'DATA_WEIGHT'},
  '1': {'data': 'DATA_ROT',
   'flag': 'FLAG',
   'id': '1',
   'uvw': 'UVW_ROT',
   'weight': 'DATA_WEIGHT'}}]
[7]:
print(mxds.xds0)
<xarray.Dataset>
Dimensions:         (baseline: 86, chan: 166, pol: 1, pol_id: 1, spw_id: 1, time: 1355, uvw_index: 3)
Coordinates:
  * baseline        (baseline) int64 0 1 2 3 4 5 6 7 ... 78 79 80 81 82 83 84 85
  * chan            (chan) float64 3.449e+11 3.449e+11 ... 3.43e+11 3.43e+11
    chan_width      (chan) float64 dask.array<chunksize=(10,), meta=np.ndarray>
    effective_bw    (chan) float64 dask.array<chunksize=(10,), meta=np.ndarray>
  * pol             (pol) int32 9
  * pol_id          (pol_id) int32 0
    resolution      (chan) float64 dask.array<chunksize=(10,), meta=np.ndarray>
  * spw_id          (spw_id) int32 0
  * time            (time) datetime64[ns] 2011-05-28T01:41:17.375999451 ... 2...
Dimensions without coordinates: uvw_index
Data variables: (12/17)
    ANTENNA1        (baseline) int32 dask.array<chunksize=(86,), meta=np.ndarray>
    ANTENNA2        (baseline) int32 dask.array<chunksize=(86,), meta=np.ndarray>
    ARRAY_ID        (time, baseline) int32 dask.array<chunksize=(1355, 86), meta=np.ndarray>
    DATA            (time, baseline, chan, pol) complex128 dask.array<chunksize=(1355, 86, 10, 1), meta=np.ndarray>
    DATA_WEIGHT     (time, baseline, chan, pol) float64 dask.array<chunksize=(1355, 86, 10, 1), meta=np.ndarray>
    EXPOSURE        (time, baseline) float64 dask.array<chunksize=(1355, 86), meta=np.ndarray>
    ...              ...
    OBSERVATION_ID  (time, baseline) int32 dask.array<chunksize=(1355, 86), meta=np.ndarray>
    PROCESSOR_ID    (time, baseline) int32 dask.array<chunksize=(1355, 86), meta=np.ndarray>
    SCAN_NUMBER     (time, baseline) int32 dask.array<chunksize=(1355, 86), meta=np.ndarray>
    STATE_ID        (time, baseline) int32 dask.array<chunksize=(1355, 86), meta=np.ndarray>
    TIME_CENTROID   (time, baseline) float64 dask.array<chunksize=(1355, 86), meta=np.ndarray>
    UVW             (time, baseline, uvw_index) float64 dask.array<chunksize=(1355, 86, 3), meta=np.ndarray>
Attributes: (12/13)
    bbc_no:           1
    corr_product:     [[0, 0]]
    data_groups:      [{'0': {'data': 'DATA', 'flag': 'FLAG', 'id': '0', 'uvw...
    freq_group:       0
    freq_group_name:
    if_conv_chain:    0
    ...               ...
    name:
    net_sideband:     1
    num_chan:         166
    num_corr:         1
    ref_frequency:    344871541760.4001
    total_bandwidth:  1864427170.9199219

Make Imaging Weights

make_imaging_weight documentation

[8]:
from ngcasa.imaging import make_imaging_weight
import matplotlib.pylab as plt

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

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

mxds2 = make_imaging_weight(mxds1, imaging_weights_parms, grid_parms, sel_parms)
######################### Start make_imaging_weights #########################
Setting data_group_in  to  {'data': 'DATA_ROT', 'flag': 'FLAG', 'id': '1', 'uvw': 'UVW_ROT', 'weight': 'DATA_WEIGHT'}
Setting default data_group_out  to  {'data': 'DATA_ROT', 'flag': 'FLAG', 'id': '1', '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 #########################

Make Gridding Convolution Functions

make_gridding_convolution_function

[9]:
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'] = 1

gcf_xds = make_gridding_convolution_function(mxds2, gcf_parms, grid_parms, sel_parms)
write_image(gcf_xds,'mosaic_gcf.gcf.zarr')
gcf_xds = xr.open_zarr('mosaic_gcf.gcf.zarr')
######################### Start make_gridding_convolution_function #########################
Setting data_group_in  to  {'data': 'DATA_ROT', 'flag': 'FLAG', 'id': '1', 'uvw': 'UVW_ROT', 'weight': 'DATA_WEIGHT', 'imaging_weight': 'DATA_WEIGHT'}
Setting default data_group_out  to  {'data': 'DATA_ROT', 'flag': 'FLAG', 'id': '2', '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  [250 250]
#########################  Created graph for make_gridding_convolution_function #########################
Time to store and execute graph  write_zarr 0.9743752479553223

Make Mosaic Primary Beam and Image

make_mosaic_pb

make_image_with_gcf

[10]:
from ngcasa.imaging import make_mosaic_pb

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

img_sel_parms = {}
img_xds= xr.Dataset()

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

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

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

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

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

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

img_xds = make_image_with_gcf(mxds2,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')
######################### Start make_mosaic_pb #########################
Setting default image_center  to  [250 250]
Setting data_group_in  to  {'data': 'DATA_ROT', 'flag': 'FLAG', 'id': '1', 'uvw': 'UVW_ROT', 'weight': 'DATA_WEIGHT', 'imaging_weight': 'DATA_WEIGHT'}
Setting default data_group_out  to  {'data': 'DATA_ROT', 'flag': 'FLAG', 'id': '2', '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_image_with_gcf #########################
Setting default image_center  to  [250 250]
Setting default single_precision  to  True
Setting default pb_limit  to  0.2
Setting data_group_in  to  {'data': 'DATA_ROT', 'flag': 'FLAG', 'id': '1', 'uvw': 'UVW_ROT', 'weight': 'DATA_WEIGHT', 'imaging_weight': 'DATA_WEIGHT'}
Setting default data_group_out  to  {'data': 'DATA_ROT', 'flag': 'FLAG', 'id': '2', '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'}
Setting default data_group_out  to  {'pb': 'PB', 'weight_pb': 'WEIGHT_PB', 'id': '1', 'weight_pb_sum_weight': 'WEIGHT_PB_SUM_WEIGHT', 'sum_weight': 'SUM_WEIGHT', 'image': 'IMAGE'}
grid sizes 500 500
#########################  Created graph for make_mosaic_with_gcf #########################
/usr/local/lib/python3.7/dist-packages/ngcasa/imaging/_imaging_utils/_normalize.py:51: RuntimeWarning: invalid value encountered in true_divide
  normalized_image = (image / sum_weights_copy) / (oversampling_correcting_func[:,:,None,None]*normalizing_image)
Time to store and execute graph  write_zarr 40.8763325214386

Compare CASA and ngCASA Primary Beams

[11]:
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
import xarray as xr

img_xds = xr.open_zarr('mosaic_img.img.zarr').isel(time=0,pol=0)
casa_img_xds = xr.open_zarr('Antennae_North_cal_lsrk_cube.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))

#print(casa_img_xds)

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)
    mosaic_pb = mosaic_pb/np.max(mosaic_pb)

    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,cmap='jet',extent=extent)
    im1 = ax0[1].imshow(casa_mosaic_pb,cmap='jet',extent=extent)
    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.m,mosaic_pb.isel(l=300),label='ngCASA PB')
    plt.plot(casa_mosaic_pb.m,casa_mosaic_pb.isel(l=300),'*',label='CASA PB',markersize=0.5)
    plt.legend()
    plt.xlabel('m')
    plt.ylabel('Amplitude')
    plt.title('Mosaic PB Cross Section')

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

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

Frequency 343.9505596639216 GHz
../_images/imaging_alma_mosaic_test_20_1.png
../_images/imaging_alma_mosaic_test_20_2.png
../_images/imaging_alma_mosaic_test_20_3.png

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
import xarray as xr

img_xds = xr.open_zarr('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 = xr.open_zarr('Antennae_North_cal_lsrk_cube.img.zarr',chunks={'l':grid_parms['image_size'][0],'m':grid_parms['image_size'][1]}).isel(time=0,pol=0)
pb_limit = 0.22
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'
#casa_image_name = 'IMAGE'

def comparison_plots(chan):
    print('Frequency',img_xds.chan.isel(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)

    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].title.set_text('ngCASA PB Corrected Dirty Image')
    ax0[1].title.set_text('CASA PB Corrected Dirty Image')
    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.imshow((100*(mosaic_img.data - casa_mosaic_img.data)/np.nanmax(casa_mosaic_img.data)),cmap='jet')
    plt.title('Difference Dirty Image')
    plt.xlabel('m'), plt.ylabel('l')
    plt.colorbar()
    plt.show()

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

Frequency 343.9505596639216 GHz
../_images/imaging_alma_mosaic_test_22_1.png
../_images/imaging_alma_mosaic_test_22_2.png
[13]:
from cngi.image.moments import moments

chan_slice = slice(0,139)
img_xds.attrs['rest_frequency'] = casa_img_xds.attrs['rest_frequency']

pb_limit = 0.25
img_xds_pb_cut = img_xds.copy(deep=True)
img_xds_pb_cut['IMAGE'] = img_xds['IMAGE'].where(img_xds['PB'].data > pb_limit,other=np.nan)
img_xds_mom8 = moments(img_xds_pb_cut.isel(chan=chan_slice), moment=8, axis='chan')
img_mom8 = img_xds_mom8.MOMENTS_MAXIMUM

casa_img_xds_pb_cut = casa_img_xds.copy(deep=True)
casa_img_xds_pb_cut['IMAGE'] = casa_img_xds['IMAGE_PBCOR'].where(img_xds['PB'].data > pb_limit,other=np.nan)
casa_img_xds_mom8 = moments(casa_img_xds_pb_cut.isel(chan=chan_slice), moment=8, axis='chan')
casa_img_mom8 = casa_img_xds_mom8.MOMENTS_MAXIMUM

plt.figure()
plt.imshow(img_mom8,extent=extent)
plt.title('ngCASA Moment 8 Image')
plt.xlabel('m'), plt.ylabel('l')
plt.colorbar()

plt.figure()
plt.imshow(casa_img_mom8,extent=extent)
plt.title('CASA Moment 8 Image')
plt.xlabel('m'), plt.ylabel('l')

plt.figure()
plt.imshow((img_mom8 - casa_img_mom8)/np.nanmax(casa_img_mom8),cmap='jet',extent=extent)
plt.title('Difference Moment 8 Image')
plt.xlabel('m'), plt.ylabel('l')
plt.colorbar()

/usr/local/lib/python3.7/dist-packages/dask/utils.py:35: RuntimeWarning: All-NaN slice encountered
  return func(*args, **kwargs)
/usr/local/lib/python3.7/dist-packages/toolz/functoolz.py:488: RuntimeWarning: All-NaN slice encountered
  ret = f(ret)
/usr/local/lib/python3.7/dist-packages/dask/utils.py:35: RuntimeWarning: All-NaN slice encountered
  return func(*args, **kwargs)
/usr/local/lib/python3.7/dist-packages/toolz/functoolz.py:488: RuntimeWarning: All-NaN slice encountered
  ret = f(ret)
/usr/local/lib/python3.7/dist-packages/dask/utils.py:35: RuntimeWarning: All-NaN slice encountered
  return func(*args, **kwargs)
/usr/local/lib/python3.7/dist-packages/toolz/functoolz.py:488: RuntimeWarning: All-NaN slice encountered
  ret = f(ret)
/usr/local/lib/python3.7/dist-packages/dask/utils.py:35: RuntimeWarning: All-NaN slice encountered
  return func(*args, **kwargs)
/usr/local/lib/python3.7/dist-packages/toolz/functoolz.py:488: RuntimeWarning: All-NaN slice encountered
  ret = f(ret)
[13]:
<matplotlib.colorbar.Colorbar at 0x7f0d0d3cc890>
../_images/imaging_alma_mosaic_test_23_2.png
../_images/imaging_alma_mosaic_test_23_3.png
../_images/imaging_alma_mosaic_test_23_4.png
[14]:
from cngi.image.moments import moments


pb_limit = 0.2
img_xds_pb_cut = img_xds.copy(deep=True)
ngcasa_pbcor = img_xds['IMAGE'].where(img_xds['PB'].data > pb_limit,other=np.nan)

casa_img_xds_pb_cut = casa_img_xds.copy(deep=True)
casa_pbcor = casa_img_xds['IMAGE_PBCOR'].where(img_xds['PB'].data > pb_limit,other=np.nan)

rmad_per_channel = np.nanmean(np.abs(ngcasa_pbcor.data-casa_pbcor.data),axis=(0,1))

print(ngcasa_pbcor)
plt.figure()
plt.plot(ngcasa_pbcor.chan[:140],rmad_per_channel[:140])
plt.title('Root Mean Average Difference')
plt.ylabel('RMAD')
plt.xlabel('Freq (Hz)')
plt.show()



<xarray.DataArray 'IMAGE' (l: 500, m: 500, chan: 166)>
dask.array<where, shape=(500, 500, 166), dtype=float64, chunksize=(500, 500, 10), chunktype=numpy.ndarray>
Coordinates:
  * chan             (chan) float64 3.449e+11 3.449e+11 ... 3.43e+11 3.43e+11
    chan_width       (chan) float64 dask.array<chunksize=(10,), meta=np.ndarray>
    declination      (l, m) float64 dask.array<chunksize=(500, 500), meta=np.ndarray>
  * l                (l) float64 0.0001576 0.0001569 ... -0.0001563 -0.0001569
  * m                (m) float64 -0.0001576 -0.0001569 ... 0.0001563 0.0001569
    pol              int32 9
    right_ascension  (l, m) float64 dask.array<chunksize=(500, 500), meta=np.ndarray>
    time             datetime64[ns] 2011-06-04T07:26:13.328236156
../_images/imaging_alma_mosaic_test_24_1.png