Source
13
13
14
14
tb = table( )
15
15
ia = image( )
16
16
im = imager( )
17
17
else:
18
18
from taskinit import *
19
19
from simutil import *
20
20
from sdimaging import sdimaging
21
21
from imregrid import imregrid
22
22
from immath import immath
23
+
from concat import concat
24
+
from feather import feather
25
+
import sdbeamutil
23
26
from casa_stack_manip import stack_frame_find
24
27
25
28
# the global tb, ia, and im tools are used
26
29
27
30
def simanalyze(
28
31
project=None,
29
32
image=None,
30
33
# if image==False:
31
34
imagename=None,
32
35
skymodel=None,
72
75
if not os.path.exists(fileroot):
73
76
msg(fileroot+" directory doesn't exist - the task expects to find results from creating the datasets there, like the skymodel.",priority="error")
74
77
# msg should raise an exception for priority=error
75
78
76
79
77
80
if not is_CASA6:
78
81
saveinputs = myf['saveinputs']
79
82
saveinputs('simanalyze',fileroot+"/"+project+".simanalyze.last")
80
83
# myparams=in_params)
81
84
else:
82
-
casalog.post("saveinputs not available in casatasks, skipping saving simanalyze inputs", priority='WARN')
85
+
casalog.post("saveinputs not available in casatasks, skipping saving simanalyze inputs", priority='INFO')
83
86
84
87
if (not image) and (not analyze):
85
88
casalog.post("No operation to be done. Exiting from task.", "WARN")
86
89
return True
87
90
88
91
grscreen = False
89
92
grfile = False
90
93
if graphics == "both":
91
94
grscreen = True
92
95
grfile = True
476
479
msg(get_taskstr('sdimaging', sdim_param), priority="info")
477
480
if not dryrun:
478
481
sdimaging(**sdim_param)
479
482
if not os.path.exists(temp_out):
480
483
raise RuntimeError("TP imaging failed.")
481
484
482
485
# Scale image by convolved beam / antenna primary beam
483
486
ia.open(temp_out)
484
487
imbeam = ia.restoringbeam()
485
488
ia.close()
486
-
beam_area_ratio = qa.getvalue(qa.convert(imbeam['major'], "arcsec")) \
487
-
* qa.getvalue(qa.convert(imbeam['minor'], "arcsec")) \
488
-
/ pb_asec**2
489
+
try:
490
+
beam_area_ratio = qa.getvalue(qa.convert(imbeam['major'], "arcsec")) \
491
+
* qa.getvalue(qa.convert(imbeam['minor'], "arcsec")) \
492
+
/ pb_asec**2
493
+
except KeyError:
494
+
# this shouldn't hit when simutil.imtclean sets restoringbeam='common'
495
+
msg("using center channel values to calculate beam area ratio",
496
+
origin='simanalyze', priority='info')
497
+
chan_index = int(imbeam['nChannels']/2)
498
+
pol_index = 0
499
+
beam_area_ratio = qa.getvalue(qa.convert(imbeam['major'], "arcsec")) \
500
+
* qa.getvalue(qa.convert(imbeam['minor'], "arcsec")) \
501
+
/ pb_asec**2
489
502
msg("Scaling TP image intensity by %f." % (beam_area_ratio),origin='simanalyze')
490
503
temp_in = temp_out
491
504
temp_out = temp_out + ".scaled"
492
505
immath(imagename=temp_in, mode='evalexpr', expr="IM0*%f" % (beam_area_ratio),
493
506
outfile=temp_out)
494
507
if not os.path.exists(temp_out):
495
508
raise RuntimeError("TP image scaling failed.")
496
509
497
510
# Regrid TP image to final resolution
498
511
msg("Regridding TP image to final resolution",origin='simanalyze')
539
552
# TODO: Define PSF of image here
540
553
# for now use default
541
554
542
555
# get image beam size form TP image
543
556
if os.path.exists(tpimage):
544
557
ia.open(tpimage)
545
558
beam = ia.restoringbeam()
546
559
ia.close()
547
560
548
561
if sd_only:
549
-
bmarea = beam['major']['value']*beam['minor']['value']*1.1331 #arcsec2
562
+
try:
563
+
bmarea = beam['major']['value']*beam['minor']['value']*1.1331 #arcsec2
564
+
except KeyError:
565
+
# this shouldn't hit when simutil.imtclean sets restoringbeam='common'
566
+
msg("using center channel values to calculate beam area",
567
+
origin='simanalyze', priority='info')
568
+
chan_index = int(beam['nChannels']/2)
569
+
pol_index = 0
570
+
bmarea = (beam['beams']['*'+str(chan_index)]['*'+str(pol_index)]['major']['value'] *
571
+
beam['beams']['*'+str(chan_index)]['*'+str(pol_index)]['minor']['value'] *
572
+
1.1331) #arcsec2
550
573
bmarea = bmarea/(cell[0]['value']*cell[1]['value']) # bm area in pix
551
574
else: del beam
552
575
#del beam
553
576
554
577
msg('generation of total power image '+tpimage+' complete.',origin='simanalyze')
555
578
# update TP ms name the for following steps
556
579
sdmsfile = tpmstoimage
557
580
sd_any = True
558
581
559
582
imagename = re.split('.image$',tpimage)[0]
581
604
sourcefieldlist = pl.where(code=='OBJ')[0]
582
605
nfld = len(sourcefieldlist)
583
606
tb.done()
584
607
elif dryrun:
585
608
nfld=1 # HACK
586
609
msfile = mstoimage[0]
587
610
588
611
# set cleanmode automatically (for interfm)
589
612
if nfld == 1:
590
613
cleanmode = "csclean"
614
+
gridder = "standard"
591
615
else:
592
616
cleanmode = "mosaic"
617
+
gridder = "mosaic"
593
618
619
+
# keep the same default minor cycle algorithm for imtclean
620
+
# as previously used in imclean
621
+
deconvolver = "clark"
594
622
595
623
596
624
597
625
# clean insists on using an existing model if its present
598
626
if os.path.exists(imagename+".image"): shutil.rmtree(imagename+".image")
599
627
if os.path.exists(imagename+".model"): shutil.rmtree(imagename+".model")
600
628
601
629
# An image in fileroot/ has priority
602
630
if len(modelimage) > 0 and os.path.exists(fileroot+"/"+modelimage):
603
631
modelimage = fileroot + "/" + modelimage
604
632
msg("Found modelimage, %s." % modelimage,origin='simanalyze')
605
633
606
634
# in simdata we use imdirection instead of model_refdir
607
635
if not myutil.isdirection(imdirection,halt=False):
608
636
imdirection=model_refdir
609
637
610
-
myutil.imclean(mstoimage,imagename,
611
-
cleanmode,cell,imsize,imdirection,
612
-
interactive,niter,threshold,weighting,
613
-
outertaper,pbcor,stokes, #sourcefieldlist=sourcefieldlist,
614
-
modelimage=modelimage,mask=mask,dryrun=dryrun)
638
+
myutil.imtclean(mstoimage,imagename,
639
+
gridder,deconvolver,
640
+
cell,imsize,imdirection,
641
+
interactive,niter,threshold,
642
+
weighting,outertaper,pbcor,stokes,
643
+
modelimage=modelimage,mask=mask,
644
+
dryrun=dryrun)
645
+
## Disabled as part of CAS-12502
646
+
#myutil.imclean(mstoimage,imagename,
647
+
# cleanmode,cell,imsize,imdirection,
648
+
# interactive,niter,threshold,weighting,
649
+
# outertaper,pbcor,stokes,
650
+
# #sourcefieldlist=sourcefieldlist,
651
+
# modelimage=modelimage,mask=mask,dryrun=dryrun)
615
652
616
653
617
654
# create imagename.flat and imagename.residual.flat:
618
655
if not dryrun:
619
656
myutil.flatimage(imagename+".image",verbose=verbose)
620
657
myutil.flatimage(imagename+".residual",verbose=verbose)
621
658
outflat_current = True
622
659
623
660
# feather
624
661
if featherimage:
660
697
cell = [qa.abs(cell[0]),qa.abs(cell[0])]
661
698
662
699
# get beam from output clean image
663
700
if verbose: msg("getting beam from "+imagename+".image",origin='simanalyze')
664
701
if os.path.exists(imagename+".image"):
665
702
ia.open(imagename+".image")
666
703
beam = ia.restoringbeam()
667
704
ia.close()
668
705
# model has units of Jy/pix - calculate beam area from clean image
669
706
# (even if we are not plotting graphics)
670
-
bmarea = beam['major']['value']*beam['minor']['value']*1.1331 #arcsec2
707
+
try:
708
+
bmarea = beam['major']['value']*beam['minor']['value']*1.1331 #arcsec2
709
+
except KeyError:
710
+
# this shouldn't hit when simutil.imtclean sets restoringbeam='common'
711
+
msg("using center channel values to calculate beam area",
712
+
origin='simanalyze', priority='info')
713
+
chan_index = int(beam['nChannels']/2)
714
+
pol_index = 0
715
+
bmarea = (beam['beams']['*'+str(chan_index)]['*'+str(pol_index)]['major']['value'] *
716
+
beam['beams']['*'+str(chan_index)]['*'+str(pol_index)]['minor']['value'] *
717
+
1.1331) #arcsec2
671
718
bmarea = bmarea/(cell[0]['value']*cell[1]['value']) # bm area in pix
672
719
msg("synthesized beam area in output pixels = %f" % bmarea,origin='simanalyze')
673
720
674
721
675
722
676
723
677
724
678
725
679
726
680
727
750
797
# get beam from output clean image
751
798
if verbose: msg("getting beam from "+imagename+".image",origin="analysis")
752
799
ia.open(imagename+".image")
753
800
beam = ia.restoringbeam()
754
801
ia.close()
755
802
# model has units of Jy/pix - calculate beam area from clean image
756
803
cell = myutil.cellsize(imagename+".image")
757
804
cell= [ qa.convert(cell[0],'arcsec'),
758
805
qa.convert(cell[1],'arcsec') ]
759
806
# (even if we are not plotting graphics)
760
-
bmarea = beam['major']['value']*beam['minor']['value']*1.1331 #arcsec2
807
+
try:
808
+
bmarea = beam['major']['value']*beam['minor']['value']*1.1331 #arcsec2
809
+
except KeyError:
810
+
# this shouldn't hit when simutil.imtclean sets restoringbeam='common'
811
+
msg("using center channel values to calculate beam area",
812
+
origin='simanalyze', priority='info')
813
+
chan_index = int(beam['nChannels']/2)
814
+
pol_index = 0
815
+
bmarea = (beam['beams']['*'+str(chan_index)]['*'+str(pol_index)]['major']['value'] *
816
+
beam['beams']['*'+str(chan_index)]['*'+str(pol_index)]['minor']['value'] *
817
+
1.1331) #arcsec2
761
818
bmarea = bmarea/(cell[0]['value']*cell[1]['value']) # bm area in pix
762
819
msg("synthesized beam area in output pixels = %f" % bmarea)
763
820
764
821
765
822
# flat output:? if the user manually cleaned, this may not exist
766
823
outflat = imagename + ".image.flat"
767
824
if (not outflat_current) or (not os.path.exists(outflat)):
768
825
# create imagename.flat and imagename.residual.flat
769
826
myutil.flatimage(imagename+".image",verbose=verbose)
770
827
if os.path.exists(imagename+".residual"):
924
981
xextent = nn[0]*cell_asec*0.5
925
982
xextent = [xextent,-xextent]
926
983
yextent = nn[1]*cell_asec*0.5
927
984
yextent = [-yextent,yextent]
928
985
flipped_array = beam_array.transpose()
929
986
ttrans_array = flipped_array.tolist()
930
987
ttrans_array.reverse()
931
988
pl.imshow(ttrans_array,interpolation='bilinear',cmap=pl.cm.jet,extent=xextent+yextent,origin="bottom")
932
989
psfim.replace(project+"/","")
933
990
pl.title(psfim,fontsize="x-small")
934
-
b = qa.convert(beam['major'],'arcsec')['value']
935
-
pl.xlim([-3*b,3*b])
936
-
pl.ylim([-3*b,3*b])
937
-
ax = pl.gca()
938
-
pl.text(0.05,0.95,"bmaj=%7.1e\nbmin=%7.1e" % (beam['major']['value'],beam['minor']['value']),transform = ax.transAxes,bbox=dict(facecolor='white', alpha=0.7),size="x-small",verticalalignment="top")
991
+
try:
992
+
b = qa.convert(beam['major'],'arcsec')['value']
993
+
pl.xlim([-3*b,3*b])
994
+
pl.ylim([-3*b,3*b])
995
+
ax = pl.gca()
996
+
pl.text(0.05,0.95,"bmaj=%7.1e\nbmin=%7.1e" % (beam['major']['value'],
997
+
beam['minor']['value']),
998
+
transform = ax.transAxes,bbox=dict(facecolor='white', alpha=0.7),size="x-small",verticalalignment="top")
999
+
except KeyError:
1000
+
# this shouldn't hit when simutil.imtclean sets restoringbeam='common'
1001
+
msg("using center channel beam values for plot configuration",
1002
+
origin='simanalyze', priority='info')
1003
+
chan_index = int(beam['nChannels']/2)
1004
+
pol_index = 0
1005
+
b = qa.convert(beam['beams']['*'+str(chan_index)]['*'+str(pol_index)]['major'],'arcsec')['value']
1006
+
pl.xlim([-3*b,3*b])
1007
+
pl.ylim([-3*b,3*b])
1008
+
ax = pl.gca()
1009
+
pl.text(0.05, 0.95, "bmaj=%7.1e\nbmin=%7.1e" % (beam['beams']['*'+str(chan_index)]['*'+str(pol_index)]['major']['value'],
1010
+
beam['beams']['*'+str(chan_index)]['*'+str(pol_index)]['minor']['value']),
1011
+
transform = ax.transAxes,bbox=dict(facecolor='white', alpha=0.7),size="x-small",verticalalignment="top")
939
1012
ia.close()
940
1013
myutil.nextfig()
941
1014
942
1015
disprange = [] # first plot will define range
943
1016
if showmodel:
944
1017
discard = myutil.statim(modelflat+".regrid",incell=cell,disprange=disprange,showstats=False)
945
1018
myutil.nextfig()
946
1019
disprange = []
947
1020
948
1021
if showconvolved:
976
1049
# if not displaying still print stats:
977
1050
# 20100505 ia.stats changed to return Jy/bm:
978
1051
msg('Simulation rms: '+str(sim_rms/bmarea)+" Jy/pix = "+
979
1052
str(sim_rms)+" Jy/bm",origin="analysis")
980
1053
msg('Simulation max: '+str(sim_max/bmarea)+" Jy/pix = "+
981
1054
str(sim_max)+" Jy/bm",origin="analysis")
982
1055
#msg('Simulation rms: '+str(sim_rms)+" Jy/pix = "+
983
1056
# str(sim_rms*bmarea)+" Jy/bm",origin="analysis")
984
1057
#msg('Simulation max: '+str(sim_max)+" Jy/pix = "+
985
1058
# str(sim_max*bmarea)+" Jy/bm",origin="analysis")
986
-
msg('Beam bmaj: '+str(beam['major']['value'])+' bmin: '+str(beam['minor']['value'])+' bpa: '+str(beam['positionangle']['value']),origin="analysis")
1059
+
try:
1060
+
msg('Beam bmaj: '+str(beam['major']['value'])+
1061
+
' bmin: '+str(beam['minor']['value'])+
1062
+
' bpa: '+str(beam['positionangle']['value']),
1063
+
origin="analysis")
1064
+
except KeyError: # per-plane beams...
1065
+
pol_index = 0
1066
+
for chan_index in range(0, beam['nChannels']):
1067
+
msg('polarization '+str(pol_index) + ' channel '+str(chan_index) + ' Beam'
1068
+
' bmaj: '+str(beam['beams']['*'+str(chan_index)]['*'+str(pol_index)]['major']['value'])+
1069
+
' bmin: '+str(beam['beams']['*'+str(chan_index)]['*'+str(pol_index)]['minor']['value'])+
1070
+
' bpa: '+str(beam['beams']['*'+str(chan_index)]['*'+str(pol_index)]['positionangle']['value']),
1071
+
origin="analysis")
1072
+
987
1073
988
1074
989
1075
990
1076
# cleanup - delete newmodel, newmodel.flat etc
991
1077
# flat kept by user request CAS-5509
992
1078
# if os.path.exists(imagename+".image.flat"):
993
1079
# shutil.rmtree(imagename+".image.flat")
994
1080
if os.path.exists(imagename+".residual.flat"):
995
1081
shutil.rmtree(imagename+".residual.flat")
996
1082
# .flux.pbcoverage is nessesary for feather.