Imaging and deconvolution demonstration

This script makes a fake data set and then deconvolves it. Finally the full and residual visibility are plotted.


In [1]:
%matplotlib inline

import os
import sys

sys.path.append(os.path.join('..', '..'))

from data_models.parameters import arl_path
results_dir = arl_path('test_results')


from matplotlib import pylab

pylab.rcParams['figure.figsize'] = (8.0, 8.0)
pylab.rcParams['image.cmap'] = 'rainbow'

import numpy

from astropy.coordinates import SkyCoord
from astropy import units as u
from astropy.wcs.utils import pixel_to_skycoord

from matplotlib import pyplot as plt

from processing_components.image.iterators import image_raster_iter

from wrappers.serial.visibility.base import create_visibility
from wrappers.serial.skycomponent.operations import create_skycomponent
from wrappers.serial.image.operations import show_image, export_image_to_fits
from wrappers.serial.image.deconvolution import deconvolve_cube, restore_cube
from wrappers.serial.visibility.iterators import vis_timeslice_iter
from wrappers.serial.simulation.configurations import create_named_configuration
from wrappers.serial.simulation.testing_support import create_test_image
from wrappers.serial.imaging.base import create_image_from_visibility
from wrappers.serial.imaging.base import advise_wide_field

from workflows.serial.imaging.imaging_serial import invert_list_serial_workflow, predict_list_serial_workflow

from data_models.polarisation import PolarisationFrame

import logging

log = logging.getLogger()
log.setLevel(logging.DEBUG)
log.addHandler(logging.StreamHandler(sys.stdout))

mpl_logger = logging.getLogger("matplotlib") 
mpl_logger.setLevel(logging.WARNING)

In [2]:
pylab.rcParams['figure.figsize'] = (12.0, 12.0)
pylab.rcParams['image.cmap'] = 'rainbow'

Construct LOW core configuration


In [3]:
lowr3 = create_named_configuration('LOWBD2', rmax=750.0)


LOWBD2
	(<Quantity -2544285.76809868 m>, <Quantity 5103070.08431238 m>, <Quantity -2848693.63167071 m>)
	GeodeticLocation(lon=<Longitude 116.4999 deg>, lat=<Latitude -26.7 deg>, height=<Quantity 300. m>)
create_configuration_from_file: Maximum radius 750.0 m includes 236 antennas/stations
[-2545983.21949368  5104856.98683873 -2852914.16896517]

In [4]:
print(lowr3.xyz)


[[-2545983.21949368  5104856.98683873 -2852914.16896517]
 [-2546185.01380468  5104746.30561217 -2852694.10384143]
 [-2546035.10056168  5104813.12493769 -2852826.9592822 ]
 [-2546265.88100668  5104811.152585   -2852823.03769538]
 [-2546083.33189368  5104781.73855973 -2852764.55441586]
 [-2546142.68682568  5104704.55486745 -2852611.09172557]
 [-2546139.85921368  5104755.56346844 -2852712.51103972]
 [-2546127.70172968  5104805.53831985 -2852811.87497164]
 [-2545996.46068868  5104756.98228663 -2852715.3320457 ]
 [-2546284.51583968  5104762.51147862 -2852726.32562029]
 [-2546228.03834568  5104761.82361752 -2852724.95796071]
 [-2545940.75549368  5104724.81905679 -2852651.38258132]
 [-2546120.75103668  5104832.97495556 -2852866.42664965]
 [-2546078.23721068  5104820.80334512 -2852842.22609603]
 [-2546168.30972068  5104779.54904667 -2852760.20105371]
 [-2546347.89420068  5104784.53197954 -2852770.10851292]
 [-2546075.02344868  5104745.83435203 -2852693.16684493]
 [-2545946.43562968  5104754.13849412 -2852709.67779366]
 [-2546171.41137368  5104809.18694873 -2852819.12946267]
 [-2546061.86944968  5104797.31197464 -2852795.5187048 ]
 [-2545935.87256968  5104778.96304978 -2852759.03592857]
 [-2545905.64079168  5104811.37478268 -2852823.47948629]
 [-2546030.10552968  5104835.44655288 -2852871.34087393]
 [-2546174.06459568  5104721.41053545 -2852644.60549099]
 [-2546000.30972768  5104719.06542816 -2852639.94276415]
 [-2546221.99697268  5104708.74471769 -2852619.4223155 ]
 [-2546107.29808268  5104727.98299268 -2852657.67336763]
 [-2546087.08128168  5104867.75586217 -2852935.58078496]
 [-2546176.84238068  5104840.10990222 -2852880.61291202]
 [-2546041.46580268  5104730.2726904  -2852662.22592479]
 [-2545995.88769168  5104824.01317535 -2852848.60813307]
 [-2546039.72190068  5104701.50119085 -2852605.02016546]
 [-2546192.67451468  5104861.94066785 -2852924.01855804]
 [-2545896.93374068  5104755.37877047 -2852712.14380869]
 [-2546214.57293468  5104819.53676339 -2852839.70777857]
 [-2546282.13665968  5104840.01027696 -2852880.41482925]
 [-2546030.49167968  5104783.05439401 -2852767.17066111]
 [-2546255.57953168  5104787.50170269 -2852776.01315016]
 [-2546037.36763368  5104869.42150387 -2852938.89254485]
 [-2546346.26240768  5104759.12982251 -2852719.60194552]
 [-2546213.60133268  5104778.85782736 -2852758.82671707]
 [-2545885.10012368  5104773.53404244 -2852748.241549  ]
 [-2546350.31075068  5104815.74541327 -2852832.16951794]
 [-2546054.53260668  5104767.41605475 -2852736.07728447]
 [-2546195.37756268  5104880.51570846 -2852960.95091538]
 [-2546316.34204068  5104738.56619846 -2852678.7157301 ]
 [-2545938.33894368  5104844.27232782 -2852888.88897413]
 [-2546319.70527968  5104798.95690768 -2852798.78929011]
 [-2546104.57012268  5104689.08779377 -2852580.33887262]
 [-2545898.90185868  5104834.14918786 -2852868.7613507 ]
 [-2546078.79774968  5104713.50937094 -2852628.895774  ]
 [-2546144.25633668  5104856.30625389 -2852912.81577285]
 [-2546256.68815468  5104736.24907566 -2852674.10864423]
 [-2545986.91905568  5104774.75361013 -2852750.66638944]
 [-2545870.36415668  5104794.07639843 -2852789.08547758]
 [-2546306.37688968  5104825.16166257 -2852850.89164573]
 [-2546070.19035068  5104845.58984121 -2852891.50855792]
 [-2546232.34827268  5104848.27412143 -2852896.84565507]
 [-2546122.93091368  5104773.02584244 -2852747.23110576]
 [-2545964.69676868  5104804.70829736 -2852810.22465561]
 [-2546277.21636668  5104861.18620374 -2852922.51847313]
 [-2546239.41107368  5104871.49628018 -2852943.01777836]
 [-2546214.39854468  5104799.76566211 -2852800.3973193 ]
 [-2546269.03248568  5104712.58140275 -2852627.05071462]
 [-2546185.06363068  5104691.84091595 -2852585.81284678]
 [-2546129.45808768  5104883.30017335 -2852966.48720756]
 [-2546300.09884468  5104782.58089536 -2852766.22921384]
 [-2545995.23995768  5104699.22380627 -2852600.49209024]
 [-2545979.49501868  5104736.27258268 -2852674.15538274]
 [-2546146.95938368  5104735.66307876 -2852672.94351909]
 [-2545942.23780968  5104825.25065334 -2852851.06858419]
 [-2545905.91483368  5104735.91494004 -2852673.44428949]
 [-2546070.88791168  5104885.34024875 -2852970.543446  ]
 [-2546328.00126768  5104840.63657193 -2852881.66007819]
 [-2546007.44726768  5104798.71512104 -2852798.30855088]
 [-2545876.77922268  5104867.81854756 -2852935.70542099]
 [-2545870.25204868  5104886.44060222 -2852972.73125535]
 [-2546354.39646268  5104667.79768964 -2852538.00821235]
 [-2545953.12473668  5104890.13568529 -2852980.07811027]
 [-2546470.39082168  5104737.56379114 -2852676.72266499]
 [-2546019.80405468  5104933.17595935 -2853065.65416985]
 [-2546284.49092668  5104895.22495983 -2852990.19700637]
 [-2545766.15358868  5104716.24682223 -2852634.33859005]
 [-2546324.52592168  5104882.96044055 -2852965.81172409]
 [-2545868.71990668  5104698.19285416 -2852598.44227014]
 [-2546225.28547268  5104907.53257509 -2853014.66797547]
 [-2546406.53911568  5104828.06646137 -2852856.66719527]
 [-2546050.72093668  5104911.17952486 -2853021.91912797]
 [-2546302.76452268  5104682.08101625 -2852566.40744622]
 [-2546482.27426468  5104758.53766905 -2852718.42457941]
 [-2546325.42278468  5104906.47028073 -2853012.55583824]
 [-2546454.84518568  5104777.30356419 -2852755.73640876]
 [-2546268.52177168  5104643.79308229 -2852490.28036319]
 [-2545769.55419568  5104849.78081115 -2852899.84137412]
 [-2546404.12256668  5104755.5198126  -2852712.42423976]
 [-2546439.86008968  5104705.5958933  -2852613.16157506]
 [-2546001.09448368  5104674.58954575 -2852551.51231509]
 [-2545788.61254768  5104761.71727629 -2852724.74652472]
 [-2545971.13674768  5104916.14230892 -2853031.78652573]
 [-2545974.67437668  5104874.62775412 -2852949.24402126]
 [-2546199.07712468  5104646.32568635 -2852495.31588586]
 [-2546247.38319468  5104676.81152076 -2852555.93022064]
 [-2545814.44720168  5104838.178407   -2852876.77256125]
 [-2545908.81718368  5104665.09493971 -2852532.63439227]
 [-2546172.71929968  5104669.67489499 -2852541.74061974]
 [-2546409.08022968  5104688.19844364 -2852578.57059672]
 [-2545815.76758468  5104738.19568282 -2852677.97904176]
 [-2546390.93119768  5104805.06817896 -2852810.94020053]
 [-2545947.34494968  5104656.44602751 -2852515.43794448]
 [-2546051.56797468  5104654.32591626 -2852511.22257248]
 [-2546346.61118768  5104857.14355221 -2852914.48055524]
 [-2546106.03998168  5104911.03792288 -2853021.63758376]
 [-2546253.97265168  5104924.54831522 -2853048.5000089 ]
 [-2546393.60933168  5104713.08232648 -2852628.0466906 ]
 [-2546465.37087668  5104845.57808792 -2852891.48518911]
 [-2546321.56128868  5104720.31130168 -2852642.41990791]
 [-2545922.35733268  5104870.46588837 -2852940.9690723 ]
 [-2546144.86670268  5104634.9012642  -2852472.60095079]
 [-2545838.47567268  5104712.75490683 -2852627.39568908]
 [-2545852.48916668  5104821.62833029 -2852843.86639647]
 [-2546409.17988068  5104731.34226058 -2852664.35252838]
 [-2546103.12517568  5104653.54066886 -2852509.6612818 ]
 [-2545897.24515168  5104713.96104212 -2852629.79382218]
 [-2545767.10027768  5104791.88352716 -2852784.72543838]
 [-2546370.49018168  5104887.99430598 -2852975.82045142]
 [-2545810.56079268  5104867.50959807 -2852935.09114329]
 [-2546422.01001268  5104886.54638404 -2852972.9415791 ]
 [-2546387.64269768  5104871.19460651 -2852942.41796703]
 [-2546458.86861568  5104825.09058165 -2852850.75031706]
 [-2546169.49308168  5104924.58133702 -2853048.56566545]
 [-2545837.14283368  5104754.6220674  -2852710.63927211]
 [-2545907.30995468  5104898.37266452 -2852996.45552052]
 [-2545823.46566368  5104687.59453689 -2852577.3698618 ]
 [-2546003.69787968  5104647.77248859 -2852498.1925319 ]
 [-2546468.60955068  5104806.24017275 -2852813.2704508 ]
 [-2545736.61937168  5104820.83356811 -2852842.28618776]
 [-2546345.22852268  5104690.6817947  -2852583.5081907 ]
 [-2546086.55811168  5104628.96909418 -2852460.80614357]
 [-2546284.36636268  5104665.01490386 -2852532.4752587 ]
 [-2546389.56098968  5104852.18356651 -2852904.6187214 ]
 [-2545823.80198768  5104887.95680626 -2852975.74589154]
 [-2546145.15320068  5104902.88321844 -2853005.42375874]
 [-2545906.69958868  5104916.2357776  -2853031.97236753]
 [-2545748.68966068  5104770.46021658 -2852742.12992655]
 [-2545844.35511168  5104852.74269817 -2852905.73043097]
 [-2545902.55159468  5104683.52893774 -2852569.28631765]
 [-2545805.87717168  5104781.51916041 -2852764.11818886]
 [-2546124.27620868  5104671.180465   -2852544.7341125 ]
 [-2545865.44386368  5104673.89888628 -2852550.13909159]
 [-2545861.42043368  5104731.62378544 -2852664.91227826]
 [-2546507.08749068  5104774.62767949 -2852750.41600424]
 [-2546519.93007968  5104830.10989453 -2852860.73010987]
 [-2545597.33147068  5104818.26178613 -2852837.17276846]
 [-2546599.43952968  5104846.64485974 -2852893.60622878]
 [-2546038.62573468  5104588.20004812 -2852379.74591823]
 [-2546010.33716168  5104971.98552526 -2853142.81840256]
 [-2546329.93201568  5104961.30157498 -2853121.57573193]
 [-2546170.42731468  5104965.52780522 -2853129.97865545]
 [-2545631.94791468  5104858.74650852 -2852917.66767911]
 [-2546548.79164868  5104880.88566425 -2852961.68649058]
 [-2546171.80997968  5104576.24559804 -2852355.97713993]
 [-2546440.43308568  5104674.03936901 -2852550.41841041]
 [-2545719.08070568  5104900.1116261  -2852999.91306076]
 [-2545857.44682968  5104953.58231008 -2853106.22768205]
 [-2545600.45803668  5104778.94457962 -2852758.99920475]
 [-2546006.65005568  5104601.36790931 -2852405.92729598]
 [-2545934.47744768  5104933.36849299 -2853066.03698038]
 [-2546538.42789168  5104730.79152444 -2852663.25751145]
 [-2545704.56895368  5104853.57272067 -2852907.380747  ]
 [-2546553.53755168  5104711.33385057 -2852624.57023322]
 [-2545901.79175168  5104945.24794454 -2853089.6566407 ]
 [-2545948.27918268  5104956.6427031  -2853112.31259628]
 [-2546105.61646268  5104998.63355535 -2853195.80211283]
 [-2545783.15662768  5104909.24970875 -2853018.08211572]
 [-2546336.06058368  5104605.48052251 -2852414.10431715]
 [-2546507.61066168  5104685.10614852 -2852572.42225223]
 [-2546245.82613968  5104969.56262109 -2853138.00099384]
 [-2546523.97842268  5104795.13309562 -2852791.18648608]
 [-2546116.07987168  5104598.0679686  -2852399.3660942 ]
 [-2545713.56250268  5104718.69155432 -2852639.19939875]
 [-2545706.99795968  5104753.97226586 -2852709.34728556]
 [-2546501.30770368  5104714.8011388  -2852631.4641685 ]
 [-2545830.86478868  5104609.89760729 -2852422.88671268]
 [-2546384.35419768  5104965.79701695 -2853130.51392339]
 [-2546439.01305168  5104905.09791718 -2853009.82719704]
 [-2546530.10699068  5104755.30321209 -2852711.99357757]
 [-2545689.70842268  5104830.8022331  -2852862.1066719 ]
 [-2545647.24442168  5104785.37375577 -2852771.78219866]
 [-2545671.32271868  5104712.3054746  -2852626.50209257]
 [-2545614.34696668  5104753.92972928 -2852709.26271098]
 [-2546066.10463868  5104955.00672499 -2853109.05981587]
 [-2546489.38689168  5104924.15485197 -2853047.7176943 ]
 [-2546211.04776268  5104586.98719685 -2852377.3344319 ]
 [-2546388.19078068  5104934.32500597 -2853067.93879475]
 [-2545858.40597568  5104624.00295191 -2852450.93206875]
 [-2545793.52038468  5104652.88471075 -2852508.35705427]
 [-2546107.47247268  5104948.50814748 -2853096.13883271]
 [-2545645.23893468  5104815.42974735 -2852831.54188613]
 [-2546541.16830768  5104666.4320425  -2852535.29292526]
 [-2545805.95191068  5104817.73063895 -2852836.11669985]
 [-2546248.94024968  5104993.53420663 -2853185.66318645]
 [-2546211.18478368  5104948.26580099 -2853095.65698035]
 [-2546631.31555668  5104784.80231052 -2852770.64600625]
 [-2546013.86233368  5104990.17549814 -2853178.98513795]
 [-2546563.83902668  5104778.83376093 -2852758.77886632]
 [-2546003.74770468  5104952.52057557 -2853104.11665797]
 [-2545930.87753668  5104574.55477021 -2852352.61530299]
 [-2545982.57175868  5104567.34930204 -2852338.2888242 ]
 [-2545924.23825568  5104601.31361945 -2852405.8193526 ]
 [-2545968.01018168  5104611.48265375 -2852426.03822677]
 [-2545966.97629768  5104587.66274437 -2852378.67760864]
 [-2545886.73191768  5104589.83882504 -2852383.00426345]
 [-2546387.75480568  5104503.70342719 -2852211.74308776]
 [-2546299.43865368  5104504.69464062 -2852213.71389631]
 [-2546372.95655668  5104480.7476818  -2852166.10066849]
 [-2546394.53110868  5104527.98508193 -2852260.02178436]
 [-2546310.08890868  5104524.48309235 -2852253.05885312]
 [-2546428.94824968  5104513.9877581  -2852232.19120371]
 [-2546643.27373868  5104801.96077189 -2852804.76180928]
 [-2546683.33364568  5104822.49417295 -2852845.58793297]
 [-2546640.18454168  5104828.8908867  -2852858.30638257]
 [-2546579.03588268  5104822.97662698 -2852846.54718603]
 [-2546688.36604768  5104848.48232687 -2852897.25962551]
 [-2546595.65277268  5104804.16819525 -2852809.1507821 ]
 [-2546586.03640168  5105023.89802846 -2853246.03492623]
 [-2546592.37672868  5105006.13006317 -2853210.70725955]
 [-2546541.39252368  5105001.18406972 -2853200.87324618]
 [-2546629.01111568  5104992.35101895 -2853183.31067962]
 [-2546644.81833668  5105019.95108412 -2853238.18730092]
 [-2546536.72135968  5105024.1247036  -2853246.48561959]
 [-2545765.68024368  5104942.40191352 -2853083.99793788]
 [-2545833.35607668  5104938.24340598 -2853075.72966597]
 [-2545723.00448368  5104964.24443289 -2853127.4269536 ]
 [-2545704.96755968  5104935.474612   -2853070.2245319 ]
 [-2545814.09842168  5104963.06852104 -2853125.08891313]
 [-2545796.18606268  5104926.4826088  -2853052.3459236 ]]

We create the visibility. This just makes the uvw, time, antenna1, antenna2, weight columns in a table


In [5]:
times = numpy.zeros([1])
frequency = numpy.array([1e8])
channel_bandwidth = numpy.array([1e6])
phasecentre = SkyCoord(ra=+15.0 * u.deg, dec=-45.0 * u.deg, frame='icrs', equinox='J2000')
vt = create_visibility(lowr3, times, frequency, channel_bandwidth=channel_bandwidth,
                       weight=1.0, phasecentre=phasecentre, polarisation_frame=PolarisationFrame('stokesI'))


create_visibility: 27730 rows, 0.003 GB
create_visibility: flagged 0/27730 visibilities below elevation limit 0.100000 (rad)
create_visibility: flagged 0/27730 visibilities below elevation limit 0.100000 (rad)

In [6]:
advice = advise_wide_field(vt, guard_band_image=3.0, delA=0.1, facets=1, wprojection_planes=1, 
                           oversampling_synthesised_beam=4.0)
cellsize = advice['cellsize']


advise_wide_field: Maximum wavelength 2.998 (meters)
advise_wide_field: Minimum wavelength 2.998 (meters)
advise_wide_field: Maximum baseline 383.0 (wavelengths)
advise_wide_field: Station/antenna diameter 35.0 (meters)
advise_wide_field: Primary beam 0.0857 (rad) 4.91 (deg)
advise_wide_field: Image field of view 0.257 (rad) 14.7 (deg)
advise_wide_field: Synthesized beam 0.00261 (rad) 0.15 (deg)
advise_wide_field: Cellsize 0.000653 (rad) 0.0374 (deg)
advice_wide_field: Npixels per side = 394
advice_wide_field: Npixels (power of 2) per side = 512
advice_wide_field: Npixels (power of 2, 3) per side = 512
advice_wide_field: W sampling for full image = 2.2 (wavelengths)
advice_wide_field: W sampling for primary beam = 19.4 (wavelengths)
advice_wide_field: Time sampling for full image = 154.8 (s)
advice_wide_field: Time sampling for primary beam = 1393.3 (s)
advice_wide_field: Frequency sampling for full image = 179176.4 (Hz)
advice_wide_field: Frequency sampling for primary beam = 1612587.4 (Hz)
advice_wide_field: Number of planes in w stack 39 (primary beam)
advice_wide_field: Number of planes in w projection 39 (primary beam)
advice_wide_field: W support = 66 (pixels) (primary beam)

Plot the synthesized uv coverage.


In [7]:
plt.clf()
plt.plot(vt.data['uvw'][:,0], vt.data['uvw'][:,1], '.', color='b')
plt.plot(-vt.data['uvw'][:,0], -vt.data['uvw'][:,1], '.', color='b')
plt.xlim([-400.0, 400.0])
plt.ylim([-400.0, 400.0])
plt.show()


Read the venerable test image, constructing an image


In [8]:
m31image = create_test_image(frequency=frequency, cellsize=cellsize)
nchan, npol, ny, nx = m31image.data.shape
m31image.wcs.wcs.crval[0] = vt.phasecentre.ra.deg
m31image.wcs.wcs.crval[1] = vt.phasecentre.dec.deg
m31image.wcs.wcs.crpix[0] = float(nx // 2)
m31image.wcs.wcs.crpix[1] = float(ny // 2)

fig=show_image(m31image)


import_image_from_fits: created >f4 image of shape (256, 256), size 0.000 (GB)
import_image_from_fits: Max, min in /Users/timcornwell/Code/algorithm-reference-library/data/models/M31.MOD = 1.006458, 0.000000
replicate_image: replicating shape (256, 256) to (1, 1, 256, 256)

In [9]:
vt = predict_list_serial_workflow([vt], [m31image], context='2d')[0]

# To check that we got the prediction right, plot the amplitude of the visibility.
uvdist=numpy.sqrt(vt.data['uvw'][:,0]**2+vt.data['uvw'][:,1]**2)
plt.clf()
plt.plot(uvdist, numpy.abs(vt.data['vis']), '.')
plt.xlabel('uvdist')
plt.ylabel('Amp Visibility')
plt.show()


shift_vis_from_image: shifting phasecentre from image phase centre <SkyCoord (ICRS): (ra, dec) in deg
    (14.94714238, -44.96258742)> to visibility phasecentre <SkyCoord (ICRS): (ra, dec) in deg
    (15., -45.)>

Make the dirty image and point spread function


In [10]:
model = create_image_from_visibility(vt, cellsize=cellsize, npixel=512)
dirty, sumwt = invert_list_serial_workflow([vt], [model], context='2d')[0]
psf, sumwt = invert_list_serial_workflow([vt], [model], context='2d', dopsf=True)[0]

show_image(dirty)
print("Max, min in dirty image = %.6f, %.6f, sumwt = %f" % (dirty.data.max(), dirty.data.min(), sumwt))

print("Max, min in PSF         = %.6f, %.6f, sumwt = %f" % (psf.data.max(), psf.data.min(), sumwt))

export_image_to_fits(dirty, '%s/imaging_dirty.fits'%(results_dir))
export_image_to_fits(psf, '%s/imaging_psf.fits'%(results_dir))


create_image_from_visibility: Parsing parameters to get definition of WCS
create_image_from_visibility: Defining single channel Image at <SkyCoord (ICRS): (ra, dec) in deg
    (15., -45.)>, starting frequency 100000000.0 Hz, and bandwidth 999999.99999 Hz
create_image_from_visibility: uvmax = 363.929962 wavelengths
create_image_from_visibility: Critical cellsize = 0.001374 radians, 0.078718 degrees
create_image_from_visibility: Cellsize          = 0.00065276 radians, 0.0374004 degrees
create_image_from_visibility: image shape is [1, 1, 512, 512]
Max, min in dirty image = 14.032041, -3.384211, sumwt = 27730.000000
Max, min in PSF         = 0.999216, -0.010932, sumwt = 27730.000000

Deconvolve using clean


In [11]:
comp, residual = deconvolve_cube(dirty, psf, niter=10000, threshold=0.001, fractional_threshold=0.001,
                                 window_shape='quarter', gain=0.7, scales=[0, 3, 10, 30])

restored = restore_cube(comp, psf, residual)

# Show the results

fig=show_image(comp)
plt.title('Solution')
fig=show_image(residual)
plt.title('Residual')
fig=show_image(restored)
plt.title('Restored')


deconvolve_cube : Cleaning inner quarter of each sky plane
deconvolve_cube : PSF support = +/- 256 pixels
deconvolve_cube : PSF shape (1, 1, 512, 512)
deconvolve_cube : Multi-scale clean of each polarisation and channel separately
deconvolve_cube : Processing pol 0, channel 0
msclean : Peak of PSF = 0.9992163257898287 at (256, 256)
msclean : Peak of Dirty = 14.032041 Jy/beam at (296, 248) 
msclean : Coupling matrix =
 [[1.         0.97145834 0.65775615 0.16083516]
 [0.97145834 0.94434792 0.6446799  0.16005303]
 [0.65775615 0.6446799  0.48814111 0.14887996]
 [0.16083516 0.16005303 0.14887996 0.09028873]]
msclean : Max abs in dirty Image = 14.043046 Jy/beam
msclean : Start of minor cycle
msclean : This minor cycle will stop at 10000 iterations or peak < 0.014043 (Jy/beam)
msclean : Timing for setup: 0.750 (s) for dirty shape (512, 512), PSF shape (512, 512) , scales [0, 3, 10, 30]
msclean : Minor cycle 0, peak [12.14574656 12.07817744 11.20193642  7.86438099] at [269, 261, 3]
msclean : Minor cycle 1000, peak [-0.23310989 -0.22818647 -0.16654131 -0.03768849] at [323, 288, 3]
msclean : Minor cycle 2000, peak [-0.00400944 -0.00213527  0.0140137   0.01272819] at [369, 189, 3]
msclean : At iteration 2001, absolute value of peak 0.012612 is below stopping threshold 0.014043
msclean : End of minor cycle
msclean : Timing for clean: 19.030 (s) for dirty shape (512, 512), PSF shape (512, 512) , scales [0, 3, 10, 30], 2002 iterations, time per clean 9.505 (ms)
restore_cube: psfwidth = Parameter('x_stddev', value=2.109653693486475)
Out[11]:
Text(0.5, 1.0, 'Restored')

Predict the visibility of the model


In [12]:
vtmodel = create_visibility(lowr3, times, frequency, channel_bandwidth=channel_bandwidth,
                            weight=1.0, phasecentre=phasecentre, 
                            polarisation_frame=PolarisationFrame('stokesI'))
vtmodel=predict_list_serial_workflow([vtmodel], [comp], context='2d')[0]


create_visibility: 27730 rows, 0.003 GB
create_visibility: flagged 0/27730 visibilities below elevation limit 0.100000 (rad)
create_visibility: flagged 0/27730 visibilities below elevation limit 0.100000 (rad)

Now we will plot the original visibility and the residual visibility.


In [13]:
uvdist=numpy.sqrt(vt.data['uvw'][:,0]**2+vt.data['uvw'][:,1]**2)
plt.clf()
plt.plot(uvdist, numpy.abs(vt.data['vis'][:]-vtmodel.data['vis'][:]), '.', color='r', 
         label='Residual')
plt.plot(uvdist, numpy.abs(vt.data['vis'][:]), '.', color='b', label='Original')

plt.xlabel('uvdist')
plt.ylabel('Amp Visibility')
plt.legend()
plt.show()



In [ ]: